{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Structured Kernel Interpollation (SKI/KISS-GP)\n",
    "\n",
    "## Overview\n",
    "\n",
    "SKI (or KISS-GP) is a great way to scale a GP up to very large datasets (100,000+ data points).\n",
    "Kernel interpolation for scalable structured Gaussian processes (KISS-GP) was introduced in this paper:\n",
    "http://proceedings.mlr.press/v37/wilson15.pdf\n",
    "\n",
    "SKI is asymptotically very fast (nearly linear), very precise (error decays cubically), and easy to use in GPyTorch!\n",
    "As you will see in this tutorial, it's really easy to apply SKI to an existing model. All you have to do is wrap your kernel module with a `GridInterpolationKernel`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "import gpytorch\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# Make plots inline\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## KISS-GP for 1D Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up training data\n",
    "\n",
    "We'll learn a simple sinusoid, but with lots of training data points. At 1000 points, this is where scalable methods start to become useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_x = torch.linspace(0, 1, 1000)\n",
    "train_y = torch.sin(train_x * (4 * math.pi) + torch.randn(train_x.size()) * 0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up the model\n",
    "\n",
    "The model should be somewhat similar to the `ExactGP` model in the [simple regression example](../01_Exact_GPs/Simple_GP_Regression.ipynb).\n",
    "\n",
    "The only difference: we're wrapping our kernel module in a `GridInterpolationKernel`. This signals to GPyTorch that you want to approximate this kernel matrix with SKI.\n",
    "\n",
    "SKI has only one hyperparameter that you need to worry about: the grid size. For 1D functions, a good starting place is to use as many grid points as training points. (Don't worry - the grid points are really cheap to use!). You can use the `gpytorch.utils.grid.choose_grid_size` helper to get a good starting point.\n",
    "\n",
    "If you want, you can also explicitly determine the grid bounds of the SKI approximation using the `grid_bounds` argument. However, it's easier if you don't use this argument - then GPyTorch automatically chooses the best bounds for you."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "class GPRegressionModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super(GPRegressionModel, self).__init__(train_x, train_y, likelihood)\n",
    "        \n",
    "        # SKI requires a grid size hyperparameter. This util can help with that. Here we are using a grid that has the same number of points as the training data (a ratio of 1.0). Performance can be sensitive to this parameter, so you may want to adjust it for your own problem on a validation set.\n",
    "        grid_size = gpytorch.utils.grid.choose_grid_size(train_x,1.0)\n",
    "        \n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = gpytorch.kernels.ScaleKernel(\n",
    "            gpytorch.kernels.GridInterpolationKernel(\n",
    "                gpytorch.kernels.RBFKernel(), grid_size=grid_size, num_dims=1\n",
    "            )\n",
    "        )\n",
    "\n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n",
    "\n",
    "    \n",
    "likelihood = gpytorch.likelihoods.GaussianLikelihood()\n",
    "model = GPRegressionModel(train_x, train_y, likelihood)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train the model hyperparameters\n",
    "\n",
    "Even with 1000 points, this model still trains fast! SKI scales (essentially) linearly with data - whereas standard GP inference scales quadratically (in GPyTorch.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is for running the notebook in our testing framework\n",
    "import os\n",
    "smoke_test = ('CI' in os.environ)\n",
    "training_iterations = 1 if smoke_test else 30\n",
    "\n",
    "\n",
    "# Find optimal model hyperparameters\n",
    "model.train()\n",
    "likelihood.train()\n",
    "\n",
    "# Use the adam optimizer\n",
    "optimizer = torch.optim.Adam([\n",
    "    {'params': model.parameters()},  # Includes GaussianLikelihood parameters\n",
    "], lr=0.1)\n",
    "\n",
    "# \"Loss\" for GPs - the marginal log likelihood\n",
    "mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
    "\n",
    "for i in range(training_iterations):\n",
    "    optimizer.zero_grad()\n",
    "    output = model(train_x)\n",
    "    loss = -mll(output, train_y)\n",
    "    loss.backward()\n",
    "    optimizer.step()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make predictions\n",
    "\n",
    "SKI is especially well-suited for predictions. It can comnpute predictive means in constant time, and with LOVE enabled (see [this notebook](./Simple_GP_Regression_With_LOVE_Fast_Variances_and_Sampling.ipynb)), predictive variances are also constant time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAADFCAYAAAB0K08/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2dd3hUVdrAfzOTMumFhCT03qXIVVAREFGkCFgogspagu4KtlUEhcVVF/jUXQTbSlFAQCMqRcEFURALEi5Ik9BLgBRSSCPJJJOZ748zk8xkZiA9E+b8nidPZuaee+e9c+97zznveYvGbDYjkUg8C219CyCRSOoeqfgSiQciFV8i8UCk4kskHohUfInEA5GKL5F4IFLxJRIPxKsqOymKEgbMBxTgdVVV42pUKolEUqtUtcdvDDwK3AmMqzlxJBJJXaCpjueeoih9AD9VVbfXmEQSiaTWqdJQH0BRlDbAPCAN2O6q3fTp06VPsERST8ybN0/j7PMqK76qqqcURbkdOKAoSqSqqmmu2v7zn/+86vEuXrxI48aNqypOneDuMrq7fOD+Mrq7fFBxGWfPnu1yW7Ws+qqqmoDfgczqHEcikdQtVbXqPwt0AX4FPlJVtaRGpZJIJLVKlRRfVdV3aloQiQTAaDSSkpKCwWCgPkLGTSYTOTk5df69laG8jBqNBl9fX6Kjo/HyqphKV3mOL5HUBikpKQQEBNC0aVM0Gqd2qVqluLgYb2/vOv/eylBeRrPZTFZWFikpKTRr1qxCx5CeexK3wmAwEBoaWi9K31DRaDSEhoZiMBgqvI9UfIlbYTabK6z0ycnJDB48mJSUlFqWyv3RaDSVmhpJxZc0WObOnctvv/3GnDlzqnyMxYsXs27dOlavXs3nn38OwGuvvcaOHTtqSkw7Ll26xNixY0vf5+XlMWHCBKZPn87KlSt58skn2bt3b618ty1yji9pcISGhlJYWFj6ftGiRSxatAi9Xk9WVlaFjxMXJ0JMRo8eDcDMmTOJiYmhVatWtWZYDAsLIzg4uPR9YGAgnTt3pn///gwYMIChQ4cyZMgQfv/9dztDXUZGBt999x0PPvhgjcghFV/S4EhISGD69Ols2LCBgoIC/Pz8GDVqFPPmzavUcb755hteeuml0vc33HADGzdupFu3bnz//ffMnz+fMWPGEBUVRXZ2Nrt27WLOnDksXryYrKwsOnbsyMGDBzGbzfj7+7N+/XpWrFhBeno6KSkpmEwmSkpK2L59O++99x4LFy6kTZs2nDx50qVMjRo1okWLFhw8eJCvv/6aLl26kJ2dTdu2bdmyZQv9+/dn4cKFtG7dmubNmzNy5Mgq/YZyqC9pcMTExBAcHIzBYECv12MwGAgODiY6OrpSxzGZTHY9u7e3NzqdDoA77riDzz77jMWLF5Obm8vmzZt57LHH2Lp1K7m5ufTq1YsLFy7QunVrbr31VqZNm8azzz7Lzp07SUtLY9SoUaxZs4bIyEjCwsL4+eef8fb2ZsSIEbRp0+aqsun1ekaPHk3jxo3Zs2cP7du3p3379rRo0YKRI0cSGhrKgQMHKvfD2SAVX9IguXjxIrGxsezYsYPY2FhSU1MrfYxhw4bx66+/lr7fu3cvI0aMKH2v1+tp0qQJN910Ew899BCxsbEUFxej0Wi46667Sttah+SjRo1i/fr1pQ+P1NRU+vXrx9/+9jdSUlLIzBQOrlcyXl66dInMzEw6dOjAihUrGDhwoJ3BMzs7m2+//ZZbbrmlWtMROdSXNEis83OABQsWVOkYEydO5P3332fdunUYjUY6d+7MjTfeiEajYevWrZw4cYJp06YRFxdHZGQkd9xxB0OGDGHRokXs2rWLZ599lhMnTpCdnc2AAQPw8fGhe/fu9OnTB4CxY8cycOBARowYwTPPPMOECRN47733SE1N5fz58zRr1oy8vDz+/PNPioqKyM7OJj4+nqVLl6LVajl8+DCffvopaWlpaLVaDh06RFZWFrt27aJDhw4cOXKEwsJC9Hp9pc+9WmG5FWH69OlmGaRTN7i7fHB1GU+cOEG7du3qUCJ7GqIDj5Xyv93s2bNdRufJob5E4oFIxZdIPBCp+BKJByIVXyLxQKTiSyQeiFR8iceyevVqunbtarce/txzz/H8889TVFRUj5LVPnIdX+KW6PW+NXaswkLn4ao333wzkZGRbN68mbvuuovc3FzOnj3LmDFj8PHxqbHvd0ek4ks8msmTJ/PRRx9x11138e233zJ06FCMRiMffvhhqT++Tqfj0qVLxMfHM2fOHB555BHGjBnD8uXL2bhxY4Wz3rgTcqgvcUsKCw019nclWrZsiZeXF8ePH6ewsBA/Pz+Ki4vt/PG7dOlCp06dSExMJDQ0lPDwcCZMmEDTpk0bbC6AKim+oihBiqKsURTllKIoH9S0UBJJXTJ58mQee+wx+vfvD0BJSYmdP/7GjRsJDw8vDQKy+s1rtVqMRmO9yV0dqtrj9wX+AnQDblcU5YYak0giqSPi4+P5448/GDx4MH379qV58+YcOHCA06dPs3XrVu6//34uXLiAwWBg6dKl5Ofns2PHDk6ePMnZs2dJTk6+YoitO1PVLLvfW18rinIIaJjjHYlHY5sJ58033yz9X94Pvl+/fnbvt23bBsDGjRtrWcLao1pWCUVRgoBEVVXPXandxYsXr3qsymROqS/cXUZ3lw+uLqPJZKK4uLiOpHGkpMT9S0S4ktFkMlVI16D6Vv2HgH9crVFFI8bcPbIM3F9Gd5cPrixjTk5OvUfH1ff3VwRnMmq12gpf/ypb9RVFGQ2sU1U1V1GUqKoeRyKR1D1Vter/DZgPbFAU5QAwvEalkkgktUpVjXsfAHIZTyJpoDQ8lyOJx/DutuovlU29re0Vty9evJiwsDD0ej1xcXF8/PHHV2xfXFzMl19+yc8//0x+fj7Lli0DRM2/p556ihkzZtCqVatqy13bSMWXeCxr1qzBaDRy//33AxAUFHTVfQ4ePEhRUREffGA/4PXy8qpw3Tp3QLrsSjyWr776qtRbD2DAgAEsWbKErVu3Mm3aNDIyMrjnnntYvXo1Q4YMwWg0sm3bNg4dOkRSUhJjxowB4LPPPmPTpk388ssvAKxdu5Yvv/ySKVOmcOTIER544AHefvttpk6dCsCKFStYvnw5c+bM4cyZMyxbtowpU6Zw/PjxOjt3qfgSj6WkpMTOZ+Ds2bOcPn2awYMHExYWRnx8vINfvqIoXHfddTRp0oSQkBAKCwv5/vvvGTZsGLfccguAXT79Zs2aERAQwAsvvMDRo0dJSEigsLCQSZMmMXHiRBYtWkRMTAyKonD69Ok6O3c51Jd4LEOGDGHr1q307NkTgHPnznH+/HkAIiIiCAkJuapffkFBARcuXADKfPit+fQ7dOjgUM7aZDJx+PDh0vclJSWEhoYyZMiQOlV82eNLPJZHH32UjIwM3nrrLVatWoVGo6Fly5asXLkSjUZD165dHfzyDx06xJ9//kl2djYnTpzAaDTSs2dPXnnlFc6ePcuxY8dK8+mvWLGCxMREkpOTOXv2LElJSTRu3Ji8vDwmTpzIqVOniI2N5bnnnuOxxx6r0xwAMq9+JXB3Gd1dPpB59WsCmVdfIpFUCan4EokHIhVf4lZoNJpaq01/LWNbWLMiSMWXuBW+vr5kZWVJ5a8EZrOZrKwsfH0rnqD0ml/OM5nMnM3M51R6PiazGa1Gg06rQe+tpUt0EMF+7m3I8TSio6NLS0rXh/KbTCa0WvfuD8vLqNFo8PX1LU0NVhGuWcXPvFzE4eRcjqbmkV/kPC+aejaLtpEBXN88hKjgypcallSNIqOJs5n5nE7PJ7+4BL2XFr23Dl8vLc3D/OrV9fVaWBmpCNec4pvNZvYkZvP76av3GGazmRMX8zhxMY8W4f4M6dIYvbeujiT1PM5k5LPj6CVyS/Iwubg26tlLRAb60rN5CO0bB6LTVnzeKqk47j2mqST5RSVsOJDCzlMZlR4mJmbm89XeJHIK6i/t07WKobiErQkX+eZAMknZBpdKbyUtz8D3CRdZtjORk2mX60hKz+KaUfzzlwr4fPd5EjPzq3yMzPwivtybxMXcK+dil1Scsxn5rN59noSU3Ervm19kZNOhFH47mYHJJI19Nck1ofiJmfms35/MZRdz+cpwucjI138kVesBIhHsOp3JhgPJ5Bmqd132JGax4WAKBUXunwizodDgFT85u5CNB1OvOnysDMUlJjYeTCU9T/b8VWXX6Uziz1yqseOdy8wnbs8FLuVf28Us64oGrfhpuQa+OZCC0WSq8WMbTSY2HUrFUCx7mcoSf+ZSjSq9ldzCYtbvT6n2CELSgK36l/KLWL8/GYPRUTFzMtNY9trTFBvySU86h7GorOfW6LwxlwgDns7bV6Qkbt6aR//5AcHhkXbHyS4oZktCGiOui6qUV5Qns/vMJXadziQnM42Vc1/gwRlv2/2uyaePsXTm4xgNhei8hcNJSbH9yMp6XRrFNMNH78ekWQtLj5FbWMyG/Snc1ysGX7kCU2WqrPiKovQHZquqensNylMhDMUlbNifQkG53vjCySO8/8KDGIsMmF2MAqxKD+KGKwEunDjM/Kfu47n3vwJg2WtPo9HApFkLOYPowfq0Dq+ls7l2OHBeLKPmZKbxn7/dy+XsTDZ9/B9uvWcS7z4/gZIiA1qdDpOlIER5hbdivS4pZ0RGmq2rP+TeKWXlGzIuG/j2UCqjukfjpWvQg9Z6o8qKr6rqDkVR/GpSmKuRnJzMQw89xCMz55Nj9rfrVQDemXKfk710wJ2A9flUZPnLA9YCIqFjXlYGr08caLen9YbbfSaLxkG+BNTGSV0jJGcX8vOJDGaM7IWxuGwevueHDez5YUPpe1MVKtXs3BjHzo1xAMxatZ3g8EiSsgrYknCRoV3laKwqVHeoXyFLS02V0Jo1axa//vorhf+Zx4jJL/Hd8nc58+ceB4UVdAQeQRT7aeLiiP8HrAP+DfzmsNV6w3l5+1C85jdub+HeM6P6KqFVUGxi/Z8ZvHzfTXZK7xo/4A5EOYYgIBO4ZPn/B7Dd5Z4bFr3J4Aen8OX8V7j/uX+hLcrn+maB1T4HK9dCGbKKUCd3cnVLaIWGhlJYWFj6Xt3yFeqWr1wcJQB4H5hk89kx4HMgB/AGfIC2wDjgXsvf78AUYI/DEae+8xmp588yYOIkftr2I927d6/Q+dQHde1uajKZWX8gGY23nhnLtrD2vdc5tPMHJy01wAOI3/wOhPK74iDwDrAaKLTbsv+nTez/aRMAv61bTkzz2VznE0TT0JobfLq7yy5UX0b37sIsJCQkMO2ll1i/fgNFhsIrtOwOfIHo7fOBVcAynPXmgmkIZf8rovL3T4gb074K6vyn7iOiaSsKLucx4cGHOXRgX3VO55pi0+4jvPxULA/OeJu5f7nTRY/fCvgYuM3ms13AeuAsEGb5awzcB1wHLAXmAf8B3gQcbTbWEdl0H1/SMzKlu3UlaBCWkZiYGPJMPldR+icQN1NHRI+hAJNxrfQAqcAsoAXiARGAuBmfcGiZfuEMACeOHUGv16PXy6CeMxn5/PvN/+PMn3vYuvpDF27STyCux23AReBviKlXX2Auold/H3gDeBrxkJgIqECkpc13gKNx1dtXT6/bhjPjk838cCSths/u2qY6RTOvA9oqitKtBuVxSmJmPv/7atUVWnwA/BfQA4uAPkCCXYuQiGh63z4SAF8/fwLDGtlszUfYA/6JMAb+F5iDGJ46EtWkGfHx8VU5lWuGkJBQOjUNZ+fGzzGbzezcGEeJ0TbOIQTYjPgtA4E4oCvwIZBc2srHz7/ckYsRD4MbgLsQD4s7EQ+CnnYtjUUG9P6BBIdHcir9MgfOZ9fkKV7TVFnxVVU9qKpqc1VVD9WkQOUpLjGx7Wg6uLTczkMM1QsQc8gnLK/tyU5PKbUuGwryybuUgUZT/vRfBR4DjMAMhNHPkdSk8/Tr16/S53Itseibn+k1cDjevs5GPsHAFoTCpgNjgfGW1/YUFeSj0Whp26MPUS3b0bJLLyKatLBs3Qz0BuKB1ojR28TSfa0PHCu/nMwgI6+I5ORkBg8eTEpKSk2c6jWJ28/xG4WFUVTkynX2JctfMWJu+J3TVr1vH8WJ/b+Tn5tNsaEQb1893W6+nfzcHE4d3E2x3RTiYyAJMeR/DmHscxxt3DBoeKXTHV0rnE6/zMUSf/7YvtHJ1mCEwt4InAIGAuecHsd6HUY8/qKdk89X7/6T9KREy7vzQH/gPeBxYCViJLYSgOhW7Xn3uQdKnXx+PJrGjmXz+O2335gzZw4LFy6sgTO+9nDrOX5aroHpn2yh18Dh6LzL5xx/AtHbm4AHKa/0Gq2WXrcNZ9aq7Yx/YQ6dbxyAsciAl49v6RDx8df/y5x1e5i1ajshEVFotFbj0P8Q800QUwdHK/6v/1uLn58foaGhNXjG7o+huIRtR9PJyUxD7xeIl49tjx+E+O36AqcR83qh9BqtFl+/ANp0v7G0te1Q3Za8rEx63z6KwFDrdMwAxAIvWt5/DAwBhJNP4pEDbPr4P8wY2YuHbmnL4sWLMZlMLFq0CL1e73HXqCK4reKbzWa2HU0nKDwCX/8ASuysxWMoq9L9JMKSX25/k8nupsrLyqTvsHFMnb+avsPGkXspo7RtcHgknW8cAGZby/FHiBvMH+HoE+bwHfqAQA7+ebha59nQWL/zMG89M5H1H86hsCAPY5F1tOSPePjeBJxBKH1i6X5mk4k3vo7HR+/PTcPH89z7XzlcByuTZi1g/47vyMsqv+1thIXfG/gKMaoQ7PlhQ9mKgmUU5ufnx/jx4zly5Ej1T/waw62G+snJyYwbNw6Ap2fNY/aMaSQe2V/OaNQB+ATxzHoJWOz0WBqN1u6mmjRrQenre6fMcmifl5Upbhg7y/RTQA/EPHMVMALbZaXCy3m0b9tGvC680orDtcGZjHwm3dXHhTv0IuAWxPLcbZb/9rw4tCte3j7M3fAH4Pw6WJmxbAsLnxlHdnpquS3TgSiEn8ZGy3ces29iNqPV6jAYDAQHB1cqF52n4FY9/ty5c4mPjyc+Pp4XpkzmzJ978PaxzRzqi+jdA4BPEU9/55jNJo7s3lHh7540awEzP/2RbjfZhh4UIpx70oGhwMsO+4U0iubHHTsr/D0NlZBQYcV3rvR/RRjd8hC/0xmnx+h0w60888G6Cn1f6SjMATNirr8JiEDYExyX+kymEkwmE4sXL5aGPie4heKHhobSokULFi1aVPpZytkTmM1mCvPzbFq+jeiBjyPWg13T67bhzFi2BYAQP296NQ+ld4tQejYPpXvTEDpEBeJVLptqcHhkuWU+EMPVByyvZwFd7LZmZ6QwZPDAip1oA+bdVetdbLkR4WUHYh6e4KIdhDVuQlBYRIW/My8rk5uGj3dyTYyI6V48Yt1/qdP9/YNCGD36nlJDn6QMtxjqJyQk8OSTT7J161ZKXAZx3IPwsitCeNflOW2l0Woxm0z4BwTRp2sbOke7ducsKCrhYFIOB85nl0b65WVlOmm5FTHnfwIxpL0V0fMIiosMhISEkJ19ba4jZ14uYtnHnzjZEgF8iXCBXohwi3bE1y8Ak6nEbuoVEehL6wh//L11FBpNFBaXYDCauHCpgFxLvL11epZ7KZ1jf/xGUYFtVqR8xDLhPmA0YtTxod335udms3bt1wAsWrSo1NjXEPzxaxu3UPyYmBgiIyOvoPQtKHuqv4gI5CjD21ePf1AI7Xr2pf89D3N421q8DTkM7nRlf2Y/Hx03tgrj+uYhHErK5deTGUyatYCczDQ+eOFhMpITbVq/BNyNmFM+SfmbbN7n2zGZzGivsaywISEhGAzOllO1CEeb5sBO4AWn+980fDy5l9KZNGsB4f4+NPUvoXf7ZgTpnd96JpOZE2mX2X8+m5QcYTeZNGsB04Zd56T1WcQoYw3CtfcXhJegc8aPH8+8efNcbvck3ELxAdLT02nSrAWFJWYyk8/ZxG1rgBUIq/oGRM9ij7HIQJc+A7l/6mxuaduIno8Nr9R3e+m09GweQqMAbzYdSiU4PBKTqfxDKBuYirAmz7PIIuqih0c3w6gP5nBKLt2aBFfuxN2cO4aP4tuvv0Cj0WK2W/X4OyLYJg0x7HbMTmwNodVoNFzfPIQ+rcPJSE9zqfQAWq2GDlGBdIgKJCW7kJ+Op3Mx18DMlT/y7eK32Lfju3J2hi8RBt5YxIhDwZkDl0ajkYY+G9xG8bdt22bXs5TFbT8KDED41T9qt89Nw8fTd9gYft+0hoKcDO6/vilRwRUvI1Se5uH+3Hd9E745kELWxWQnLb5GLO3dg3AouQeAzJTzvDi0K2g0HD9xkuZNXYUBNxzKR0TaK30n4DXL60lYH4B2aDQEh0cS7u/D4M6RVSpYEh2i575eTdh+PJ0EwNc/AMzmUl+MMp5BjMS6APMRIzJHUlPLrxB4Lm5h3AN4c+Umeg0s31M3Bt6yvH4WsF/X3b3la5q06cTj015n28Z11VJ6KxGBvoy5vik6L1fPxCmI3n80wuJvg9nMU8+9VG0Z3IGEhAQGDB3txCVXi1hO1SP8HBy9JXU6Lzr27ken6CDGKU2rVaXIS6dlcKfGDOwQyWXLkqvRwZOzAOESXIiwwwx1OI7ZbObTVZ9VWY5rDbdR/EKfUPFEt2M+Yoj/P8objqxW+yC9N6N6xODnU3MhmYF6Lw78eZjIpi2cbE1CrCWDsGbbB5ls2bDm2vAW8w8jr8hUzp0ZhBtzX4Qr7fNOd9V6efH24lUM7hRZY6mxrmsazDdrv+SNz7a7iBE4SNly64fgJF9ScKB/g78uiecvMGbMmGovT7qN4gPs+u5Lm3dDgAkI6+1fHdrq/QOJjo5hdI8YAq8wZ6wqbVo0w9fls+QjhA9/c8rcSO1pyNVeTSYz24+l26XMEnREhM+CCHm2XcXQ0LZHH24aPp7eNw9kcKfIGo9jiA7R8+BtPfALCHTS64Ow/6hASxs57SksLKR///4Ndl3/mRmziY+Pr/bypFspvlZn1TQ/ylxyX8WZQ8jOjXE8P6w7of61V+22d69e3D/xEWLnLLGJGAOxlPes5fU0wLHI447dB2pNrtogOTmZ/v37079/f0JCQ3jk1nblWmgRQ3s9Yqhffohv5sl5HzP99f/jx41ray14KSrYl9+/W+PiwVqCMPIZEbEWNzg9Rk0oTl0TGhqKXq/nuzUrMZvN1Y5DcAvFT05O5r8vxzL1Hetw/h9AG2A/YrjviFar5ejRo7UqV1xcHCuXfsioYXfSrmffclt/QcSY+yNy95Wh8/KuUNFOd8LWa7Jn/6FOwqCfAm5GGPKcD/FfHnU9d3RqXOtLmidPnmToqPtKh/tePr6ENIqyhFnvQyztaYEluLJfN7QAnt1/HKT3oBGl51zdOAS3UPy5c+dy+s+9liy57RE3lglhqHFePGHChAl1tjQzsH0ExXlZKEPuI3bOEpst0xCGpQkIpRCUGIv5fNE7HEqqfL24usbak9h6Te7eur5czEIMZUPnvwGODjAajZbDCQl14scQExNDs8bhdtGW+bnZNisPryKyJ3fHlX+BlYYQwGMymdl/SQsaLcWGQnx8fKodh1Cviu/sphPJL3wQw8ldTvfr3Lkzubl1p1ReOi2b1n/F6Cem88lsW1fhRIQbMcACbDP27NwYxw1tGxMS4t49SkJCAiNHjkSnu5JxdD4izn4dwn/BkfvHjq/TZcyLFy8yefJk3lm5gagW7SgusjVCFlC2pDcbkVjVOa1atXL7Xl9NzOLo6XPs3/E/AO644w5iY2OrtTxZr+v4CQkJtG3bFlOpQ8adCO+4HOAVp/totTo6dOhAXFyc0+21RYifNz2bBjqJGpuHSNulAA8Dy0v30Wg0LPn2lzqVs7LExMQQFRV1Ba/JOxEu0pcpy1FgT8s2bSkqrNty1nFxccLX4KOPXLTYinD8ehix+nK301bt2rVj69attSNkDZCUVcig7q3sjJkbN4oEKNXJ+1ivPX5MTAwPPGANgPGibD7/BsJhx56R947h1KmTda70VrpF+9O2RbNyUWP5lC3vzcV2GSkgJJzjFy+TlOXeIbuu6x7oEYkwQQyfHTPptGzbgZ7XdauXa5KQkMC4cePw8XXlv/EiYuVhhOXPkRMnTrhtr19QVEKHFlEuVjCqN02pTrLNvyuK8pCiKFOq/O1AXl4e7du3RwzNugAnEMNmezQaDdERYfXqcqnVaBjUKbIsdr+U1YhpSQzCp1+Ql5XB/Kn3s35nglvXd5//0XJ63z7KyZbpQDvK8tzbo9Fq6dG1c709iGNiYggODsZY7OguLLiIMBSDuKecPyDcMVmH2Wwm7udDxLTuSEflVqdtqvPAqpLiK4rSD2ikquqnQJiiKH2q9O2IIdvx4xmIDLcgfMAdc7ObzWaWLFni8Hld0zjIl3cWLWfWym00irEu8ZkRji0gjEnNS9vnZqbx3D03ERIa4pZx4SEhIbSJDmXPD+XDbttTNpL5K+WNrBqNlj8OHeWLL+pH6a1cvHiR2NhYIiIiCAgKoVmHbuUeyu8jHlxtcOVz4R8QiNlsZvDgwRWq+lQX7D6bxfIP/sP5Ywc5c/gPh+06na5aD6yqzvGHURZ4fdjy3rklDtdDyfbt21v8899FJFPYiivj0ejRo5k5c2a9XhhrOGebADP7g4PLRe/tRHgXjkcM+R+027e4qIhffvmFmTNn1toacmXCTVNTU5kyZQrX3zyAndu2oNHqMNsFJn2A6CGXAr867H/b0LsJD/St9PWo6ZDYd999F4BXXnkFg9HE8h1HeOMvd9m0KEG4Wf+E8Oz7FNvsQI1iWnDg2OnS8mxvvvkmb7/9NvVJu3bt7RLMGuxyUjQD0jGZDHh5eaHVaqukE1VV/AhEsTMQDtJXHH+7Kvdz5MgRWre+EZFRpYQypxhH1q1bx//+9796j6W2nssInyAWKbeSkXSW9KRziF7/JYQP/0SEF5lj7v2VK1eycuXKWosLr2hppdatW9sZ9OyVfjwwGBEb4Tz2INCr6mWcarNE1YkfxMjR1z/QRmF2IKZjExBr/GXFVTOSE8lITkT9Wbz/4osv+OKLL9/uoZIAABmCSURBVOotbj/PYOTVlVtZ88E8J1mM/bBWeTKb72blypVVHgVXdY6fRpmTehDlo2cqSExMDEOHKoicdk8Bfzptp9PpaNq0qVvNw1o28udfH35qceyxzt8TETcWuHI88tWXOV7UR/536xKqayt+MGXnMI3yl9bHz5+Q0FAWLnS0w9Qn1vP6/FORMMS+lwQxBctFBFbdSenSq6W2gtaSYVmv19fbnN9YYuK7Q6l4B4Wzb4ezVPELEb4JvkAmPj7lM09XnKoq/ibKck53QUTRVAkfHx+8vI4h/N9tPtf7o9Fo0Ov1mM1mhg8f7nax1Le0DRcRY3bMBVIQDj3j7LZ4eYsyYF56f6Kjo5k7d26dp4X66aefCA8Pv4JL7esII+WvCF8Ke4oK8snNyXE7l1erhd/PT2Rb0mi1+PrbVtFNpiyU+F1Epl5KMytb8y8UFhbWW9z+tmPppclHzA7G4AcRI+MCRP6DPIqKiurWuKeq6q9AoaIojwBZqqpWPKtlOeLi4izOI/Y3olarYfLkyezYsaPazgq1RZi/D/MXLefZ974irLHVeSUPmGl5/Ra20XtT3/mMvsPGsWr5slLHpbrO/7506VIyM125E/dCjLyMCIOe85UId8xZb7XwW3M6mE0mJ73+AoRpqgPl3Y41Wh2dbujP8FH31su9tjcxiyMpuSL704uTbOJWADojSpGBSAYjsgxptdo6N+6hqqrz8KcqkJ2dTVBwKMVFhbTvdRM5mWlgyGPBAjGctP53R25oGcqRTl3xtnOm+AShOL0RCTpn0LZHH5q06cTgCU9y/sRhmjVtwr6dP1FQUICfnx+jRo2q1bRQ5RNrOKJF3Fw6hPek8xRWfn5+dSZzZbFa+EeNGsXUqVM5efJkuRbFCEPfD4iH8yqsvglmUwlH1V94dP1uHh/YqU7lPpuRz2+nxMhx6+oPOX1IFQVkSkB0HGsoyyxdlli0Om7r9Z6Bp/wNefwPkaq6oVSj9fXW0bdNOAV5uTZ5+U0In/adiOXJFZzcv4sXh3blpuHjOHf0ACmnj2EsNqDX6+sk/3tCQgJTn36ajRs3ukiR/QQiY+55hLOOPd4+vkQ0Cic1NbXOZK4stv4EgwYN4uSpU2jQlMse9CMiRftYxANuLAAhETFEt2pHflEJmw+nMqp7TJ3EHWReLuKLnw/xjwcG2I3CygrIvIcoNppA+fD06rit13uQjnVuZo068vbVM2bsOLcy5F2NrjFBLFy/k1krtxESEWX5NB6RC86bMu83Sos8FhcVivThhYVMmDCh1oeXUVHRFHkHu1D6FpRFGD6NswzGxUUGfH19iY2Ndevpl5WLFy/y4MSJBDudivwd4YI8BhB1FLLTkzmx73cAzl8qYOdpZ9mWa5bsgmLW709m3eK3XUy9HrH85VtkLXOL1mq11XKcqnfFL/W+som0CgsNcaue5GpoNBpubdeotAhEmeHsZcQCyG2U5eZ3ZPXq1bXq/WY2m/n+yEWSUlwp6keIxZkvETkFnXPmzBk++ugj+vfvz4IFC+rNY68ixMXF4e/vT/alTGJatUPnZZu34TzCiAmiRxXbSozFvHr/jeRkprE3MYuTabUXf5BXaKRZVAR/HdTRScITELZza4fxN8qveHl7Vy8PRb0rPoinc5+77mf6+1/wuJv3JK5oFuZHm4gA8rIyuX7QSMunmYglMRBLZM4z8JpMplozlJnNZn44msax1DwmzVrArFXb6XbT7TYlwh9G1KHPQMx/BX6BwXhZqhhZI/caSi268lGfyWdOlCvDBmK59Sgiceg0uy2vP3gbAFsTLpJ52dGLtLrkF5Ww8MsfKHZZBToIMa/3Q8zpl9tt9ff3r3YuCrdQ/Li4OO55cjpj77yFdxcudOue5Erc0rYRj/xjIT529onliKWxaMp6GUe69hnA7n2Halym7cfSSUgumwtaqwWJeW80ZT74z2IbGFWQl1MaHGJd8y8oKHC7eb0zyi/t+fn50W/IKEIibOUuomzOPAthObdgNvPi0K78fXgPVvy4n4GDbq8xX4vC4hLW7Uvmo9fEqoLzpK5LECsP+7F9GFvJz8+v9jWod+OelchAbzpEBV69oRvTqkmkE8u5GXGD7UVcxLXAdod905PPs+2cEe+gfFo18nfYXllKTGZ+PpHBoaQcu89njOxVVlWWDxDJTDdirTfvDJ1Ox+DBg4mMjGwQozHbpT2rIbJt0wiyc7txXqcjLyvDkkR0G2U5+ZcgKiSV2UDMZjPrl73Hrp2/MfPV11jy3w+cfl9FuZhroEVkiN1nJcbyiWamIAyOOcD9CMdYe2oirZlb9PgASrOGrfQgepoxY8a6yAD7BuLnXgE4DulTE0/y9B2d6NYqptppuzLyiliz5wIHLziW9Joy/zP8g0IRzkX3IG4wxzz0ev9AOweqli1bsmTJkgYzGrMu7VkNkVkZ6SxYvJyXl22hY+9+NtGILyIyJ9+M8GEoo8RYzM6NcZjNZlYu+7ha07Ht+44zcNDtxP5riY0BuDy3IFYaQNSQOOHQIqJJC+LjHV3BK4vb9Pihfm4jSpWJiYkhNDQEY5GzeeEbiLl0X0T6Z0djn0ajYcayzew+c4nUHAO3d4ysVAZhs9nM3sQsfj+dSYmTMOCczDQWvxJLfm4jRA1AEK6s5x3aFubn8cQTT/DYY4+xdOlSt4sqvBq2DyirH0iR0cTR1DwS4n+ymfNnI0Zk6xFelxuwBvFYbRzGIgPevnquu3kw/5o7t1Kl0opLTPx0PJ1/vf4Gpw/t4dBv31NUmO+kZXNEwRYfhP3hK6fHM5tMREW5enBUnIavbW6G6GkeJ6jjzSz/9z9IT7JG8JUg3C73IYJgvkU4kJRhNpuZ+5c7mbvhDxIz81m4cQ9fvv0iS5etoFfHVi6/s7C4hDMZ+fx25BKXzc6tvWVDfD0iCjIYsZ692KFtkxatuK5zxwbhQFUZfLy09GvXyMnwegNlkZWLEKndsUuAYSwy4OsfwKFLWs7Fn+P6FqF0igp0WTcgt9DIoaQc7ujZ2u441uVce/wRD57GwBachQ+HRkbTvmdfGvk4z0FZWaTi1zDWnuZsRj7bN31to/ggEkA+jUhT/T4iU681RFRDUHgEYY1jyMlMIzg8ku9Xf0DCvt08//JsYl/6Fy0b+eHrpcVHp8XHS0tSUjLP/+1xxk57i6CwCC5fLiIgQCh+TmYay157Go0Gkk4esZnXv4NwzT2OmNs6MuSOwXz4/ns1+8O4AVd2CnsaUQvwTpxV3gVNabXf7IJith1NI/70Jbo1DcbfR4dWIxK1mM1wMv0yZzLyMZvNzPhkM98ufotDO3+g2FAoqjmbzeWSmX5C2TUZj8Vlzx4zlGSn8t8vaqYakFT8WqJlI3/OJTgmUBAXeQQiSmwVMAhhYTbjHxjMuaMHeeOhQXaONjs3xrFzYxxe3j7M3VB2zK/f+xdH9+/m+1UfcO+Uf9h9y9bVH3LuqMjtHxAajjErExGW+gTCYDQGMb+3Z/i948hMT6vGmbsvu3btYsyYMSQmJjrZmoZYL49DPBz3APFotFp6DhjKiMdfJDg80m6P5JRk3npmIg/OeNthm5Xg8Eh8/QMwFhnQefvYeORZeZkyY95IyqLdBQ/O+DcnD+zm1L5f2afuYs6cOcycOZPqIhW/Fri6X3wswj32FoSx7wHATGqixbe8nGHP21dPt5tvZ8TjYghob5m3fzAAdtsASwRhJ8oiIJ9GLBXZE944hi9XLau1Yhj1TY8ePQgIcCytVcYXiGvyNMKZqTdmUxr7tn/HhGlvkpOZxsq5L5Qq+tbVH3Lmzz1sXf2hw4PXlksXkwkMi6BZ+64k7Npus+Vh4F+IlYQJgKN/xMq5f7d7v2jRotIAqerkC3Abq/61REJCAlrtlX7aTGA4wrA0jrL4d0e8vH0wFhnQ+weW9iozlm2xqx/n7aun123DeeaDdcxYtoWOSr9yR2kJbAYCEct2jvN6H70/N97Q+5pVeitZWVl06dKFZSs+pXGz1qUGvDJeQPhdNEfM+3WYgXemjmHdh3M5fUjl9YkDeXFo11KL/86Ncbw4tCszRvZy+p3H9v5KbmZaOaUfg5jygXAhLp90Q4QWA7Rp09bOJ2H06NHVdqKSil8L2GcPdsUBxHJaEcJ55u9OWxmLizCbzez/ebOIWsR++Gh1c9b7BxIUFgHAUdU2pXczxHp1C4RN4Qmn31NUmM/2H9w3zXRNcfr0afbu3cv4sWP4/tfdKIPFsl6ZJ2MxQilTENOwN8Bs4sKJwxz8ZbPTY1ofvDOWbbH7fMbIXrw4tKuT+IgRiGmeDpEM1DGRKVC636lTJykoKADAYDAQGBhYbQceqfi1RF5eHp07d0aj0VyhF92GGO6BKMwxweXx8nOy2Lq6zOCUl5VJ32HjmDp/NX2HjSMzNYlXx/Th9YkDbfaKQUSjtUakRByGCPhwpDqx3Q2VjlGBlFzOthjcbJUzGTESMyISjj5yxeMUGwrZt/07h3l+lz4DnbQejJhGeCNqMth7c2q0Otr26ENYVNPSUaNWq6Vdu3Zs2rSJ2NhY0tPTK3GWzpGKX0vExcXRoUMHJk+ezL333nullpRl6P0UYexx/qCwDimn392TyzlZDJ7wJE3adOLeKbM4tvfXcraBxoi48/YIr8G7EKmnnFOXJcnciQ1r19B7kLNiGzsoG4V9zNVKcZktSUsunDzCi0O78uLQrhz4ZUu5VuMRy3a+iDRaMxyPYyrh1IHd3NhvIFC2EjFo0CAGDRrEggULylWeqhrSuFeLbN68+SpGPivvIAIy3kAYe/ogRgLlPO80Gjr2voULJw6XGpV2b1nrYMwTSr4UaILwGrwTp/XutFrMJhNdunSp05Jk7sLVjbALEYa3BYhsSlGIgB5742tYVFOCwhqRk5nGZ29OcziK6N3/jcieA8LI6jqxbNcb++FTnEtsbGytOVBJxa9FEhISmD59Ohs2bKCgoAC93g/fwGByMtOczPvmAn8g5n4jEXXe70PYAiyYzXbzd0dnkADElMHqgvsrYtnQSS5UjYZXPlpLSvxG0i6mNhhX3JrEen2ufO7vIZb6ViB6/SjgMYQtQHAp9QKXUi+Um2ZZaYaItOsLGBAK/18n7co4vm8XanbZQ782HKjkUL8WKR8sUlRk4Obb7mDmpz/Sa+BwJxbl/yHSde1FVLDZi3DjvO0q3xSEGEbuRyi9AdEz9UdUk3HkH6u288jd/Xnv3YYbDVldrNdHo9FcZRUmDrEKkwc8hPChn4oYpbmiOWL+vg+h9GeBflxN6aNimtR6+XeoYo+vKIofYmJqUlXVfZKuuSHWYBHrkC0pOZlmTZrg6x/gxJkD4AxlwRqPIyz/9yDqlqxG9D45iGlAU8u22ykrD7UPcXNeOcT3tQkDeLOecse7E9aqu2lpaXz99ddoNOUNfVa2AgMRPX8XxDRgluX/IYQh0IgYdT2MeFBYE2ZuQlyT8ll9NFinDdZp16i7R9SJraVKiq+qaoGiKCq2ReElTnEWLJKcXci7Fqt80umjnHUokVSIiBR7DZhMWV1BV/lNTQhj1BrE/NFVLTmBVqtl7NixbpUos76wXp9x48YRGzuZjz/5hBKjM8UH4c3XDTEVexnhhOUqx0IRYqTwAc4qEWm0OgJDw2nZqQdBYRHcN+Fhjv20rs7Cnqszx69wapKKlPhpCD1PTcmoA954eyH7kkRuu28WzWPv9yLllX04birixpqLqNDTBxFcE2L5X4hw/NiAGAlc4Tst7qKiRzNXq/xSdXDX62wtxTV09P28Ou/fHPzlexc9vxlhmV+PmIL9BTHV8kaokwaRb2Eprq5Jq24Kf3m1LLY/1M+LEV3C8bm1G3B1famJ3/Cqiq8oysuIdCC2rMOZmdgFFS2ZVJullWqKmpLxzohIDLoUEjPzMeTl0HfYOPoOG8Pvm9bw+3dflDP+GRFrv19W+fuCwyPpfMOtzHz+KdbHrSQlJaXefm93vs7KddChRQwHMZc6R7lmm+Wv4vj6+fPUW2WptHy9dIzt3ZRQ/8rl0Kvub3hVxVdV1WnJFEVRBlbrmz0crVbD0K5RrN2XxKRZZVbbe6fM4tDOH7iclYnWy4sSY7GLzLhOsKT31mi16HTelJQY0el0BIdH0qRtZ955ZwHdm4Uw6OYbaumsrg2K8y4xfOxDtO8/iveen0hJcTHtevbh2N7fqnFUDRqtFh+/slgBH52WUT2iK630NYFczqtHfLy03N09hjV7LpBTWDYv/8eq7VU63uXLl10GofRpHU73ZiFOt0nsiYuLo8Rk5us/kpizbk+tfIe3Tlz7qOD6qR9RpeU8RVG8EIa9roqihNWsSJ6Fv4+OUT1i8PPWXb1xFVFahnFjK3mZKoPOMiLz96n5vlGn1TC8WxRNQuuvaExVrfpGwL2qJjZgQv29ubt7NOv3p2AwuqpiWzV6NQ/lpjbhNXpMTyFQ78Xd3aPZsD+ZguKauS5ajXigNA+vfkLVaslRr98uKSUqWM/Y3k1pFFDeqafq9GweSr92jWrseJ5I4yBf7ru+CUH66s/D/X28GNUjhtYRV8oJUDdIxXcjQv29GdO7CR2jgqp1HL23jmHdorlVKn2NEObvw/3XN6nWQ7lJqB/jlaY0C7uSt1/dIRXfzfDWabmzS2MGtI9AV4Wija0jAph4Y3PaRtZ/r3ItEejrxX29YqqkuNe3COWeHjEE+LqPLd19JJHY0b1ZCO0aB3I4OYdDF3LINbjOrqpBQ1SwL02jvLi5i+eF1tYVvt467unZhLMZ+aiJWSRlFbhsq9VoaNc4gO5NQ4gJcb/Kz1Lx3Rh/Hx1KyzCubx7KqfTLnM0swGA0UVRioqjYhI+3ljYR/rRpFECg3qvOvfA8lZaN/GnZyJ/k7EL2nc/mssGI2QwmsxkNGtpEBtAlJgh/n9pbqakuUvEbAFqthnaNA2nXuOFXG7qWiAnRu2VvXhHkHF8i8UCk4kskHohUfInEA5GKL5F4IFLxJRIPRCq+ROKBSMWXSDwQqfgSiQciFV8i8UCk4kskHohUfInEA5GKL5F4IFLxJRIPpNLReYqi3A68iqjf9JSqqt/VtFASiaR2qUqPH6yq6q1ALKI0q0QiaWBUWvFVVV1rebkbSK5ZcSQSSV1wxaG+q/JZqqquA4YBFaq6KGvn1Q3uLh+4v4zuLh/UQe28K5TPigACVFX9vCJfImvn1R3uLh+4v4zuLh9UX8ZKD/UVRQkAhqmqulRRFC9FUWQOZ4mkgVEpq76iKL6IusxBiqI8jajXfH1tCCaRSGqPSim+qqoGYGDtiCKRSOoK6cAjkXggUvElEg9EKr5E4oFIxZdIPBCp+BKJByIVXyLxQKTiSyQeiFR8icQDkYovkXggUvElEg9EKr5E4oFIxZdIPBCp+BKJByIVXyLxQKTiSyQeiFR8icQDkYovkXggUvElEg9EKr5E4oFUuoQWgKIoQ4ApQDQwQlXV1BqVSiKR1CpV7fETVVW9G1gL9KlBeSQSSR1QJcVXVTXB8jIX2Fpz4kgkkrrgqkN9V2W0LPu+DJyzvHfJ7NmzqyqfRCKpBTRms7nKOyuKciMwW1XV4TUnkkQiqW2qa9U/BfxZE4JIJJK6o9JWfUVRNMAm4AcgG3i9poWSSCS1S7WG+hKJpGEiHXgkEg9EKr5E4oFUyXOvuiiK8nfgIhCiqup7Np93AMYB+cA3qqoeqw/5LLK4kvEB4FkgGHhIVVXVneSz2b4U+FRV1e11LZuNDC5lVBSlE3ArcEhV1Z3uJJ+iKPcAjSxv81VVXV0f8llk6Y9YObu93Oc3A7cgOu9PVFW9WJnj1nmPryhKP6CRqqqfAmGKoth6/i0A5gPvAfPqWjYrrmS0GDbzVVXtA7wN/NOd5LPZfjcQWB+y2cjgUkZFUToCsaqqLq5Hpb/Sb/iMqqpLVFVdAjxaH/JZUVV1B+DnZNNcxD34GVW4D+tjqD8MsHr+Hba8R1EUP6Ctqqp5qqoagNaKotTLiMSVjKqqmlVVXW/5fDeQXA+ygQv5ABRFaY0YySU42a8ucSkjsBA4qyjKAosC1gdXkm+PoiivKYqiAB/UuWSOFNm+sYyMjZb7MRExcqoU9aH4EcAly+tCRKAPQBiQY9POCETWoVy2uJLRlsHAf+pMInucymd5UA5VVXVtPclliysZA4BWiFHdv4E1iqL4uIt8FmYBbYG3gB11LFdFsJUdhO5UivpQ/DTA3/I6CMiwvM4A9Dbt/IGsOpTLFlcyAqAoSjvgrKqqh+taMAuu5OsPPKgoynbgL8A7iqI0rXPpBK5k9AEKVFU1WXqrJJw/WOtLPoB/AU8ihtKf17FcFcFWdgBDZQ9QH4q/Cehued0F2KwoSohleH9WURR/RVH0wDlVVQvqQT6XMgIoihIF9FBV9StFUQItPZhbyKeq6o+qqt6squpAYBnwrKqqF+pBvivJeAkwKIpitUGkAfUho8trDHRXVTVXVdWNgHc9yOYURVF0iqIEqap6HEsnqShKG2B7ZY9V54qvquqvQKGiKI8gevQs4L+WzS8B04DngOfrWjYrrmRUFKURsBmYoSiKCvyEWIFwC/nqWo4rcRUZpwCzLSsk/6eqaombyTdfUZSpiqLcC3xU17LZoijKdUBbRVG6IewQMy2b3lAUZRrwIPBKZY8rPfckEg9EOvBIJB6IVHyJxAORii+ReCBS8SUSD0QqvkTigUjFl0g8EKn4EokHIhVfIvFA/h94sKW4I4rWqwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Put model & likelihood into eval mode\n",
    "model.eval()\n",
    "likelihood.eval()\n",
    "\n",
    "# Initalize plot\n",
    "f, ax = plt.subplots(1, 1, figsize=(4, 3))\n",
    "\n",
    "# The gpytorch.settings.fast_pred_var flag activates LOVE (for fast variances)\n",
    "# See https://arxiv.org/abs/1803.06058\n",
    "with torch.no_grad(), gpytorch.settings.fast_pred_var():\n",
    "    test_x = torch.linspace(0, 1, 51)\n",
    "    prediction = likelihood(model(test_x))\n",
    "    mean = prediction.mean\n",
    "    # Get lower and upper predictive bounds\n",
    "    lower, upper = prediction.confidence_region()\n",
    "\n",
    "# Plot the training data as black stars\n",
    "def ax_plot():\n",
    "    if smoke_test: return  # this is for running the notebook in our testing framework\n",
    "    \n",
    "    ax.plot(train_x.detach().numpy(), train_y.detach().numpy(), 'k*')\n",
    "    # Plot predictive means as blue line\n",
    "    ax.plot(test_x.detach().numpy(), mean.detach().numpy(), 'b')\n",
    "    # Plot confidence bounds as lightly shaded region\n",
    "    ax.fill_between(test_x.detach().numpy(), lower.detach().numpy(), upper.detach().numpy(), alpha=0.5)\n",
    "    ax.set_ylim([-3, 3])\n",
    "    ax.legend(['Observed Data', 'Mean', 'Confidence'])\n",
    "\n",
    "ax_plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## KISS-GP for 2D-4D Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For 2-4D functions, SKI (or KISS-GP) can work very well out-of-the-box on larger datasets (100,000+ data points).\n",
    "Kernel interpolation for scalable structured Gaussian processes (KISS-GP) was introduced in this paper:\n",
    "http://proceedings.mlr.press/v37/wilson15.pdf\n",
    "\n",
    "One thing to watch out for with multidimensional SKI - you can't use as fine-grain of a grid. If you have a high dimensional problem, you may want to try one of the other scalable regression methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up train data\n",
    "\n",
    "Here we're learning a simple sin function - but in 2 dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We make an nxn grid of training points spaced every 1/(n-1) on [0,1]x[0,1]\n",
    "n = 40\n",
    "train_x = torch.zeros(pow(n, 2), 2)\n",
    "for i in range(n):\n",
    "    for j in range(n):\n",
    "        train_x[i * n + j][0] = float(i) / (n-1)\n",
    "        train_x[i * n + j][1] = float(j) / (n-1)\n",
    "# True function is sin( 2*pi*(x0+x1))\n",
    "train_y = torch.sin((train_x[:, 0] + train_x[:, 1]) * (2 * math.pi)) + torch.randn_like(train_x[:, 0]).mul(0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The model\n",
    "\n",
    "As with the 1D case, applying SKI to a multidimensional kernel is as simple as wrapping that kernel with a `GridInterpolationKernel`. You'll want to be sure to set `num_dims` though!\n",
    "\n",
    "SKI has only one hyperparameter that you need to worry about: the grid size. For 1D functions, a good starting place is to use as many grid points as training points. (Don't worry - the grid points are really cheap to use!). You can use the `gpytorch.utils.grid.choose_grid_size` helper to get a good starting point.\n",
    "\n",
    "If you want, you can also explicitly determine the grid bounds of the SKI approximation using the `grid_bounds` argument. However, it's easier if you don't use this argument - then GPyTorch automatically chooses the best bounds for you."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "class GPRegressionModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super(GPRegressionModel, self).__init__(train_x, train_y, likelihood)\n",
    "        \n",
    "        # SKI requires a grid size hyperparameter. This util can help with that\n",
    "        grid_size = gpytorch.utils.grid.choose_grid_size(train_x)\n",
    "        \n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = gpytorch.kernels.ScaleKernel(\n",
    "            gpytorch.kernels.GridInterpolationKernel(\n",
    "                gpytorch.kernels.RBFKernel(), grid_size=grid_size, num_dims=2\n",
    "            )\n",
    "        )\n",
    "\n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n",
    "\n",
    "    \n",
    "likelihood = gpytorch.likelihoods.GaussianLikelihood()\n",
    "model = GPRegressionModel(train_x, train_y, likelihood)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train the model hyperparameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 57.8 s, sys: 847 ms, total: 58.7 s\n",
      "Wall time: 10.5 s\n"
     ]
    }
   ],
   "source": [
    "# Find optimal model hyperparameters\n",
    "model.train()\n",
    "likelihood.train()\n",
    "\n",
    "# Use the adam optimizer\n",
    "optimizer = torch.optim.Adam([\n",
    "    {'params': model.parameters()},  # Includes GaussianLikelihood parameters\n",
    "], lr=0.1)\n",
    "\n",
    "# \"Loss\" for GPs - the marginal log likelihood\n",
    "mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
    "\n",
    "def train():\n",
    "    for i in range(training_iterations):\n",
    "        optimizer.zero_grad()\n",
    "        output = model(train_x)\n",
    "        loss = -mll(output, train_y)\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "\n",
    "%time train()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make predictions with the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAO0AAADPCAYAAAAZKVayAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAZ2ElEQVR4nO2deZRcZZXAfzchpNNZm4ROlGGyuoEgwnUiCoKCy+B4BnHJcIRBBBVRRHFkkRGOOjhxRRRROYoLjuMRMcg24DaIuHKNyOoICURUQkhMSDqQpLvqmz/e61AU/b5X9aq7ql7n/s55J6/e97333X6pW/db7nevhBBwHKc8TOi0AI7jNIcrreOUDFdaxykZrrSOUzJcaR2nZLjSOk7JcKV1nJIx7pVWRF4lIlek5+8Skde28KzFInJbi/KcICIXtfKMuuftJiJ7R8oPFJFPiMgiEVnZjHwi8h4ReU16fHsUZP2giJwpIvNFZNx/98aKrn9xInKIiGwVkfNE5DsiclaTj7gJ2DM9vySEsCKjnRkiMiv2oBDCKmBb3X37i8ifRGT58BdRRF4sIjeLSP8Ij/kVMLXJv2FE0vbeAqwVkTNE5CERObhO5pXA2SGE1cBjDTy2Vr7PhhCuAX4CzBkFkX8H7A78BThdRCaOwjN3ObpeaUMItwCPAB8HjgfOE5HnNnH/tprzaqTqh4Go0qY8SWlDCLcDnwZ6ap4fSBRl3Qj3Vxpoo1FOAX4ZQhgErgZWhRB+WV8phDDc5lADz9wpX83fs71VQVO2pc8dIvkheMcoPXeXouuVto4pwCCwRUT+S0Q+IiK3icgUEXm7iLxXRL4iCUtF5N0i8u8AIjI9veew9PMpIvIWEblKRPYCFHi9iMwSkeNE5M1p2RwR6ReRc0Xk7cCzR5Dra+m9venng0MIv0jb+ayIHC8iX621LCKyUERWp+eHi8iP0vNXiMibROTb6d/wHBF5p4h8UUSOrGv32BDCHbEXJiIfqO+diMgnReRt6Xs7JpXveyKyuKbOjFSGF6eXekTkfSLy2+EfzVSuE0TkUhFZGLn2FhE5Djh5+PkhhN8Dx8Zkd0amTEr7GuA04A0hhDXAauA24CCSL8ME4B5gB0lXbjnwORIrSAhhC3A/ICJyFLAphHAZcD2wAbgP+C7w98CLgLXp8/cBPgJ8N4TwJeBP9YKFEDYB1wFvEpEZwKaa4vXAt4B5QH/NPffXPOsWEsEmAqem8twM7AccDPQAp5N0K0nrTgH2aOC9rQQm1dx3FHBVCOFSYAbwOpKezK/S9obl2ww8AEh6afcQwqeAbwCvFJF9geeFEL4OXA5cmHFtIXBoCOGb6XuoZXb6vpwm2K3TAjTBNSGEK2o+V4ANIYSKiOwDfCuE8DPghnQs2ReS3RCPiQx/73Z2D/cjUXDSLy81dZ4DPBRCuCF91gTgEmDYWmWNCz9P8oUeBGrHzbcBbyL5Mckbw80DpqRtD49ZpwJfBn4OvLGm7hxgIOd58NQu8VkkPwAAS4CBur91Uca9m9N/twJzgQNq2v89sG/GtaUkP0Lw1Hf3KMkPx2achimTpY1xH3CGiEwUkZeQjMGeLiKz0/L6v/Ne4F9h5+zq00jGoZI+68S0WzyXxJL/lcT6kt7zlPeWjm03A/uHEP6W1usDzgshfIPkCyt1tw2JSA8wO5VxHaAicpCI7Aa8AnhWCGEZcCnwrpp71wMzs15IOkQ4ZISiU4Avp5b6fuC1IrJARKYBL8l6Xv3jgbtIhhSQ/LCszLj2V5LewjC1724m8HCDbTopXa+0IrKUxKq8ruZaL4m1PFISE/klkvHufcCSEMKjJBblahE5FZghIvNJfvlfAHwfGBSRu0i6bg+RWMQPknRZrwTuAM4NIdwKnA18RETOJelqHpYh7sXAVTWfHyMZC342/Xw0yZf6WalCf4eke/4aki/wPBLFvI7EWv8aOExEPgH0kXTfAQghPA48kP5Q7Qa8Eni2iJyVjmGvJ/0RAJ4rybLQImAyiRX8Jokl/VjazhdIrPmwfHuRDA1eALwQ2Cvt6u6Xvsc/AD8RkfOBNwNnhRBuG+HazcDdIvJd4FXAIhGZmv7w3ZdOojlNIL6ftryIyD8Bj4cQftxpWZol7RHNCCFc22lZykbXW1onm/QLPyu12qUh7Yrv4QpbDLe04wARmZHO9paCssnbbbjSOk7J8O6x45SMUVunPfvss91kO13D8uXL65fXsutecFrYtCXqDr5m+fLlC1oWapQYVeeKb++d7VO+b89U7tq2NbP83S+7seX239P3QLR87cZ9mNd3d2b5XTseb1mGk+4+Plr+jMos7p24KbN86Ht7ZpY1Qv/374uWz99/T9bc/ki0TmX9+pZkqB5yQLR84aJp3L867hey+rWTC7e/bG3876tn05apfOj072SWn3/RG+cXFmYMKJNHlOOMGVXK01F0pXUcYDCM5uarscWV1nEYh5ZWVd9H4hI308wuHluRHKf9DBLbat1d5C75qOohwGwzuxzoU9WlYy+W47SXSgiZR7fRyDrtUaTb2IC708+OM64YJGQe3UYj3eM5wMb0fBvJTpQR2bcne63r73ePT+H3DixoQJQ4a+mNlm8ayIx/BsDWodajqjyjEo9Y87TqtGh5pT/+N+Qxc//4klH/gszdfDupbi6+3AIQFsX/xrnzpuQ+Y8qUSbl1RpPB7tPNTBpR2kdgpzZM54kNzU8htg6bV37EtAcaECXOvJx12qRO9jrthlFYp7134vMbqBNZp13X2pe1P2cNFhj7ddqe/LXm3HXa5xX/4Wg4gFgNladsde5eGukeXw/sn57vA9wwduI4TmcYDJJ5dBu5SmtmPwe2qeqJwCYzu3nsxXKc9lJBMo9uo6ElHzP7j7EWxHE6yWAoz94Zd65wHKBSog1vrrSOwy5saWM7dXoHFkRniN85a1XL7d+zI75ks3VoOxt3ZGfGeNsfjmtZhqEV8ZnTSn9vdIZ47rWrW2t/w9+i5dUtPVRy6oQXPa8lGVYfE5/5nTJlEqsPiNc56Yj/Ldz+wDebnz+u7KpK6zhlZTA3JHX34ErrOMBgKKa0WX75qnojSQqZAAQzW5hev4okJO01ZvbWIm2Wp0/gOGNIhQmZRxZZfvmqOh04w8zmkyjuf6fXXwB8wczmFVVYcEvrOAAMhkKqMJJf/q/NbAtJtgVIskT8ID1/KXCaqv4EeIeZNZJ69Cm4pXUcoBIk84jQiF/+ocDPAMzs48BCkpQuZxeV1ZXWcUgsbdYRIeqXr6q7ARUz2xkWw8yGSFLWLCwqqyut41BsTMtT/fJvVNXabVQvBXauXanqsNmeTpretAiutI5DMnucdWRR75efHl+sqfJSkoz3w9yiqhcDx5CkLy2ET0Q5DsWXfEbwyz+2puwDdXVfXKiROlxpHQf3iHKc0lHU0nYCV1rHwXf5OE7pcEvrOCWj6mNaxykXbmkdp2TsskobSzW5lt5oiNO8DeyNcHLOJvYllT7ueyh7g/f2K+e2LEP/1fHN/DP325P+O7JDmA6tay5NYz15G9jDoumESbOjdVa9vqclGU55+Q+j5dMG5nPotDXROmf03Vu4/Q8WCKLqSz6OUzJ2WUvrOGWlGtnN02022JXWcYhb2taSpIw+rrSOQ9zSdhuutI5D3NLmpwtrL1GlTWPdXAYcBNxgZqe2RSrHaTND1dEN7JaWPSmIm6o+E1gGPJZe+2ORNvPG2C8E3kySiOyINDCV44w7qkjmkUUs4XpGELeLgAuBi4HlRWWNWloz27ngpqp3Amtj9ddu3CezrB25YZdU+qLleblhB1vMDQvJOmyM/gUzouXVR3dvqf2waHq0vJHcsD0t5oadNjA//vxt/bnPeJg256ctZmlHDOyWfn5SEDeSUKqLzWwAQFUXqupuafiZpmhoTJt2k/9kZg/G6sVyv+aVxyL/N0rMcWJnnYkbM8u2r2tNYYCo48QwayJ1Kq06V+w+J7fOA6u3RMtXHdiac0We4wTAQE6duS04V8C+Td9RcCIqM7CbmX1cVT8NfIwkiNsXgc019w4BewIPNdtoo0tQxwPnNftwxykLQ2FC5hEhGtitLojbBqD217CXJDxN0+QqraoeDVxlZltUtXU/P8fpQqphQuYRITOwW30QNzPbDqxR1V5V7QEeNLPHi8galUhVTyUZOF+tqrcDry7SiON0O0UsbU5gt5GCuJ0FnAm8FzijqKx5E1GXAJcUfbjjlIWhajFnxazAbiMFcTOzO4E7CzVUgztXOA7uEeU4pSNnwqmrGFWlvWtH9rh669B2NkTKRyOhc95+2MH+3uiyTt5e2EbIW7KpPrp7tE6rCZ3z9sL2TJmUu6STtx82j/f0xR191jGJ/pw6q4a2tSRDs7ildZySUXRM2wlcaR0Ht7SOUzo83IzjlIyKd48dp1x499hxSkbc0lbbJkcjuNI6DhBCpyVoHFdaxyFvIsotreN0HT6mdZySUa0WU9qsGFGqeizwHmAGcLyZWXr9c8AbgNvM7FVF2izPPLfjjCGV6oTMI4usGFHpXtrHzGwp8EngQ+n1vYCVadyoQgoLrrSOAyQTUVlHhJFiRGFmwcy+n16/lSdCyrwM+KCqXquq+XGBMnCldRygWp2QeUTIjBFVw5HApwFSi7wY+PHwtSK40joOSajErCNCNEaUqi4B1pjZzoiGqRW+ECgcRdCV1nGAUJXMI0IsRtRc4HlmdqWqTlPVqcNxo1R1d5JucyFcaR0HCEEyjyyyYkSp6mzgRuAcVTXgpyRZBb6jql8DTuKJWFJNM6pLPifdfXxm2TMqs7h34vMzy4dWxIN8N0KrCZ1bjTkMDSZ1jsQmbjWh8ztf/oNo+dSB+RyeE3P49L77WpLhD4M7ouVbK4Nsyqnz9j9kf5fyeEmBe4ou+WTFiAIOGKH6Gwo1Uoev0zoO5HWDuwpXWseB3BmnbsKV1nFwS+s4pSM24dRtuNI6DkCJlLbhJR9VfbaqXjeWwjhOxyjoXdEJGk11ORl4BTB1bMVxnA4xDse0J5IkETomVukZlVmZZXkJnSvjIKEztJ7UudWEzlNzEjpPbiihc2vvYWtlMFq+fWtcRoDFOQnC4zSfQTJ01z73KLlKq6pHAj8zs8dUNVr33onxlxUrH1rXeubvTid0htaTOrea0DnPcQJga25C59acKzYObs+tM3VmPA/VqgYShGexFwWs5jgb074V+Lyq3gQcoKrnjq1IjtN+pJp9dBu5ltbMlg2fq+pNZnbB2IrkOB2gRJbWl3wcB7otdluUppTWzA4fIzkcp7N04dJOFm5pHQeQ0Q/s9kxgGcmWvGvM7I8jXSvSpu+ndRwo5FyRFdgt5SLgQuBiYHnkWtOMqqUdujJ7nbTS3xtd1pl73erW2+9wQmdoPalzqwmdz9gj/h7XSg/z+uJ17tnRWkLnk++J74VdUunjvr+OtN30CXasyF9PzmTWuqZvkWLd45ECu/1aVacAi81sAEBVF6rq9BGu7WZmQ8026pbWcSDxiMo6sskK7NYHbK6pN0QS/7j+WqHIDz6mdRwoOhGVFdhtA1DbneoFBka41rzrFm5pHQco7FwxYmA3M9sOrFHVXlXtAR40s0dHuPZ4EVldaR0HCk1EZQV2S4vPAs4E3gucEbnWNN49dhyKL/lkBXYzszuBO+vqPuVaEVxpHYfu9DHOwpXWccA9ohynbLildZyy4ZbWccpFQY+ojuBK6zjgltZxyoaPaR2nZLjSOk7Z8O6x45SLXdbS9l+dHXpz5v570n979l7WoQ1/a7n9TueGhdbzw+bth83jnh2PRcu3Dm1nY06dt+Tsh81j8Mr4XthKfy+DOSFz514TzzUc5YR47OkRcUvrOOVil7W0jlNaYpa2ib0EqtoPvAtYC9xmZr+oKZsIfAk4lGTjwLFmtiONavE7YBpwjpl9NdaGb81zHEY1WPlHgW+a2SXAOapaq/JLgXNJ9t7uAfxzev144B/MbF6ewoIrreMAiUdU1tEkrwDurfm8YPjEzH5hZg+bWYXEsj6UFh0G3J1GdszFu8eOA4WClavqB4Bn1l3e08yGVX04btT9dfdNBKaa2S2QZPFQ1T5ghar+zsx+Emu30VSXApxAEt/192b2l0buc5yyUMT32Mw+Wn9NVV9c87E2blQt/wJ8uO5ZG1X1TOAQIKq0jXaPlwO3mtn1rrDOeGQUu8c3qeqS9HxyGqR8empdUdWXAivN7C+qOje9Njzu7QNuymugkVSXB5MMoP+sqscB55vZjqb/FMfpZkZvyed84DRVXZueA5xHoswTgc8Bj6TnK1T168C1qvo9EsO4Mq+BRrrHRwOXmdk3VPVLJNPZnx6p4vz9s8O49i+YGW2kuqV1x4ZOJ3SG1pM6r5XW3sPWoXhu2G0NJHRe0lJC5/wE4U+fOTn3GXkJwuM0H2x9tLbmmdlfgXPqrr2/5uPVI9y2XzNtNKK0PTwRZPla4LVZFddEPJ7yyiuj4RE1aXZunbFM6AytJ3XOi/6fR563E8DUmXdEy/Oi/+eR5+0E8Md1cTkbSRCeRd+BzXtElcm5opEx7S3A89PzScCtYyeO43SIAiFUO0Wu0prZFcBUVV0GzAe+MuZSOU6bGVeZ4AHM7N/GWhDH6STdqJxZuHOF40BXdoOzcKV1HECq5dFaV1rHYReOxlhZvz6zrLp5crS8GxI6521gb4RWkzrnbWLPI28DeyMJnfM2secx99r4O5i13xzm3pH9XQAYerj5xNBPML6XfNzSOg74mNZxyoZbWscpGT4R5TglY5ediHKcsiKVTkvQOK60jgOjNhEVC+yWli8AfkHiQnycmf1IVU8GKsAc4FNmFh1he4woxyEZ02YdTRIL7AawDJifBnH7UarEL0kDuj0MvCGvAVdax6E9gd1Udfe0fI2qvmmE+ncBR+U14N1jx6HYkk+zgd3SiC9HqOrfAdep6q0kXeKNdfWjuNI6DsWWfIoGdjOzP6vqBcBzgUdIYkNl1q/Hu8eOA6O5CT4vsNvwGHcK8CvgRmDf9No+wA15DbildRxAKqO2UBsL7PY34POqegXw8zSeFKp6q6qeRNI1Xp7XgCut48CoLfk0ENjtwBHuubiZNlxpHQd3Y3Sc0rHLujFWD8nepxkWTaPakx3LdvUx+bFw8zjl5T+Mlk8bmM+hkfClp/Xdm1nWKHftiMdx3zq0nQ07Hs8sP7kbEjrn7IfNo7JhY7S8uqUnt05r+6sHmr7DLa3jlIxRnIgac1xpHQd8E7zjlA3vHjtO2QiutI5TKsbVmFZVe0kWi1eSpLz8qJltjt/lOCWjPDrbkO/xK4H1ZrYCeBA4YmxFcpz2I9Vq5tFtNNI9/g3wIVW9jmQXwo1ZFRcumpb5kLzcsFNGITfstJzcsD05uWHX0boMWyuD0fLtOflhuyE37Kz95rQkQ16u4f7Fe+Q/Y3HzsYufoMA67XjqHqdp5i8CLgUuN7PMaNr3r46/rFj56gNad66IOU4MMxCp09/3x5Zl2DQYd64AmDrzzsyybsgNmxdIPI88xwmAB377ULS8OmFW4fbnPK3ATeNpIkpV9wb2Av4RuFFVV5vZT8dcMsdpJ6NkaWMxolR1PvA7kiTtk4Cvmdm5I8WNirXRyJj2IGCjmW0HPsMTCaYdZ9wgIWQeTRKLETUbmGdmC4CzgO+l158UNyqvgUaU9gZgb1U9CngW8PUm/gDHKQfVavbRHJkxosxsZRpyBuBAM/ttRtyoKI2MabcBZ6Yfr2/koY5TOgqMaZuNEVVzXy+wFUaOG2Vm0ckVd65wHIrNHheNEUUyP/Q/dc+qjRsVVVqPEeU4AJVq9tEc0RhRKS8Efjn8YYS4UVFcaR0Hku5x1tEc5wMnqerpPDlG1KtgZ+zjHcNdaFU9GPitqp4D3D8cNyrGqHaPYxvZp0yZFF2LfduRP265/ffvsSpavlYmM68vu86qwezN6Y3yjv+Lb2JfPNTHqrX7Z5bvWNFiQudr4u9g5n570n/HI9E6rSV0zt/AXl08PXcd9v6j404iMXR9884V8QmniZGyJ5MXIyodw55b8/mXjBA3KoaPaR0HwLfmOU7JqMbS5nWXmnSXNI7TKdzSOk7J6MLdPFm40joOQKU8WaVdaR0HxtcuH8fZJWjeiaJjuNI6DhCCK63jlAu3tI5TMnz22HHKRfDZY8cpGd49dpyS4RNRjlMuvHvsOCUj7Kq+x8seiu/TfG6k7NHLs/eYNsoHaOQZz2m5nRiH5NZ4lKfFYg/Mam0vKyfkBfneRt+BeXVaCRQO+cHCB3JjExfaE1ucNQe9f2Esinx+QO02IqFE7luO43i4GccpHa60jlMyXGkdp2S40jpOyXCldZyS0ZZ1WlV9H7AOmGlmF7ejzbr2pwOXkSQTu8HMTm23DKkczwY+ZWav7lD7ApxA8n/xezP7S5vb7yUJL7oSWAp81Mw2t1OG8cCYW1pVPQSYbWaXA32qunSs2xyBFwJvJlkqPkJVX9BuAVR1MkmipantbruG5cCtZnZ9uxU25ZXAejNbATwIHNEBGUpPO7rHRwH3pOd3p5/bipn90My2pgmx7yTJHdpuTgS+3IF2gZ2R7JcCL1PV/0wj3beb35BE319C4sFxYwdkKD3tUNo5wHBq8OEsYh0h7Sb/ycwebHO7RwI/S380OsXRwGVm9jlgD5LEx20lte4XAZcCD3f4fZSWdijtI8BwjoesLGLt4niSvCrt5q3A51X1JuAAVT03p/5Y0EOSgRzgWuJepWOCqu4N7EWSNe4EVT2s3TKMB9qhtNfDTqfgfUiSVLcdVT0auMrMtqjq3Ha2bWbLzOxwMzscuM3MLmhn+ym3AM9PzycBt3ZAhoOAjWa2HfhMjTxOE4y50prZz4FtqnoisMnMbh7rNutR1VOBC4GrVfV2oCOzt53EzK4ApqrqMmA+8JUOiHEDsLeqHgU8C/h6B2QoPb5hwHFKhjtXOE7JcKV1nJLhSus4JcOV1nFKhiut45QMV1rHKRmutI5TMlxpHadk/D8HV2pxzGwsGgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAO0AAADPCAYAAAAZKVayAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAW1UlEQVR4nO3de5RdZXnH8e+TQO4BwiUEARMaikK482hAsIAoWrpWRbswxUKRImpRikK5V8W20GgLigaKFBFElxWsoAJCUUTlIvIURAErYGu4mHCxgRAgIZl5+8feQw7D7L3PZebM2TO/z1p7MWdfzn4T5sn77sv7PJZSQkTqY8JoN0BEWqOgFakZBa1IzShoRWpGQStSMwpakZpR0IrUjIJ2CGZ2ipmd0cHxM83s62a27zC26TVmNqVk+3lmtpuZHWVm57fSPjP7vplNNrMlZvbeDtu5lZn90My2MbM/6OS7ZGi1D1oz+7yZWcU+E8xsmxa+9r+ASQ3Hm5mdYGZrzWzXfN1EM7vIzE4ffHBK6Tngt0Bpu5plZjsCOwCb5QH2tSH+PCenlO4F7gSml33fEO17e0ppDfArOvydSCktA36ff5xoZn/WyffJq9U6aM1sHrAIeFvFrn8DbN/CV/c1fkiZ84H7yf/OUkp9wMPA4ma+o0MnpJRuSSk9DvwU+GFK6bFBbRw437omv/Pl9jUcu6bjlmZW59/7ELCHmc0epu8Vah60wDuB44APNa40s780syPM7Hoz2wjYH3iHme1nZr/M9znCzC7Le+FL8s8XVpzvgvx8mNl0YGVKKZnZXDO72Mw+amaD2/L+/PunmtmF+XkmmNmxZvYhM/uGmU3K2/xeM7tx0PF7AS+WNcrMZg8MSRvWzTCz75nZPmY2zcw+YGYnmdnFg47d3cx+3LBqgZldYGbX5ds3MbOz8uPPzUcdQ62bbmZ/Z2bvA/Zp+L4Ajqr4e5UW1DZozWwyWW9xDbCbmW2Vr18AbJVS+ipwOTANuBe4IaV0K+uHbrfm/51KNiy8Eti74rRfAw4xs02AdwNX5+u3BL4FXAW8a9AxtwCklF4E7s7XvROYQzZEfQzYLl93D3DKoOPfCCwra1RK6UngiYZVE8n+IXtPSukO4IP5+vuBfjPbtOHYn/PK34NHU0ofBl5jZpsBZwA3pZQuBjbJ2znUuhOBO1JKl5H9fQ94glcGsXSotkELHEp23Xkk2S/JMfn6ncl7ppTSN1JKy8u+JKX0PNkv1qFUDA9TSi+QBffRwBYppafyTfcCO+bLxCbavhPwcErphpTSScBDwLlk/wj81aBr9M2B55r4zsZh8Y7AnwIvNJzvl/n5PgQ8U3Lsyvy/zwNTgD0azn8vsKBg3b6s/4dj4Lzk59qoifZLk+octNumlM7L/2U/nuyXfSJZABxuZhua2WvNbBcgkd1PmgBskB+/GTDBzHYCFqWUrgTWVd3UIhsin8ore5MPAv3AT4bY/yWy3v7lc5JdC3/EzKbkI4O5+X67ArsBuzQc/wSwKQXMbEFjz5m7D/ghcGb++WHgxHxYfkBDe6oYWe/s+efpZKOFodb9DnhTw7EDv1ub5ttkmGxQvUvvMbMPANub2ZSU0mpgQ2AG8A/AWcAdZL9Y30wpnWFm25MF9iPATWZ2PXAj2RD1KWAnMzuHrPd4N9lwdRczm573xC9LKf3GzK4Cbm5Y/d/AErLh+kb5PwQLgGeB24FNzOzLwAqyX+JPkPXsDwFLUkqfNrN/A74O/Az4dcN330w29CS/BNgTeGM+RJ9B1ot+EJgPvAHYOP/5LOAXZrYMuIhs6P4w8Cmym1kLgGfNbDkw38zmkwXiC2Z2G9mQ34Gzgc+b2QbAupTS98wshlj3a+Dfzey1wBZkQ+JHgG2AH5T875QWmebT9j4z+2fg9JRSs3eGe4aZnQV8Jr+0kGGgoK2B/IbQH6WUrq7cuYfkI44pKaW7K3eWpiloa8LMJpH9/xquZ6kjzsw2SimtrN5TWqGgFamZOt89FhmXhu3u8WmnnaYuW3rG4sWLm37ve/HZx6dnnit9XXvp4sWL53XcqGEyrI98vjFni8JtC6ZO4/4Xi28gfvjg/+z4/Cdu+j+l25ev2Ik5sx4o3P7g2ucLtzXr6F8dWbp9+75ZPDxxReH2Nd/csqPzz/72w6Xb5+66BUt/8VTpPn1PlW+vkt60W+n2efNn8tvflL8v8pvDprZ9/kXLWmv/M89N51MnXFm4/ZPnv2du240ZAbV8Tisy3Pqpz0BRQSsCrE3DOSlrZCloRRiDPa27nwQ8CWwcEUtGtkki3beW/tFuQtMqH/m4+37AZhFxBTDL3ReOfLNEuqsvpcKl1zTznPYQsvmmAA/kn0XGlLWkwqXXNDM83pxsdgpkaUTmFO24YGrxjK/XTp5cepLpqzq/q768OO8ZAM+s2rZ0+6p1qztuw/Z9s0q3b9U/o3T72tnNzpob2sa7Fj92A5g9b+PK7+h/dlLlPmXS/Jml27ecU/04Z8rUztrQqrW9F5uFmgnap1g//3Im6zM/vErZc9iq7QfMWNpEU8rNmVX+nDbbp/g57cpheE778LLdq/cpe077ZGe/rLMrnsECI/+cdtLmlftUPqfds/3ntDu3cUzf8OTg64pmhsfXk03Ohmzu5g0j1xyR0bE2WeHSayqDNiJuA1a7+9HAMxHx46pjROqmDytcek1Tj3wi4h9HuiEio2ltqs/cGb1cIQL01WjCm4JWhHHc05bN1Jm+am7pHeKqGTrNqJqls2rd6tI7xFUzdJpRNUtn7exppXeIq2bpVKm689v/7KTKfapm6VSpmqEzZeqkyrvDH35b+7O+nv7KHi0f0zdeg1akrtY2la66NyhoRYC1qb2gLXov391vBF5PlnM7RcR2+fpryCpZfDcijm3nnPUZE4iMoD4mFC5Fit7Ld/eZwIkRMZcscL+er38D8K8RMafdgAX1tCIArE1thcJQ7+XfGRHPkSXLBzgYGLhAPxA43t1vBv46ItrKBa2eVgToS1a4lGjmvfw3k5eLiYjPkFWveBo4rd22KmhFyHraoqVE6Xv57r4B0BcRL6fFiIh1ZLWgtmu3rQpaEdq7puXV7+Xf6O6N06gOJCuEBoC7D3TbM1lfarVlCloRsrvHRUuRwe/l58tFDbscyCsLtd3q7kvIirxd0m5bdSNKhPYf+QzxXv7hDdvOGLTvvm2dZBAFrQh6I0qkdtrtaUeDglYEzfIRqR31tCI1069rWpF6UU8rUjPjNmjLJrIvtymlKU67VmayJMVpp2UmoXoS+8a7blGa5nSky0ym+TMrU5x2UmYSqiewVyVEgM6SIpyBJsGLjHnjtqcVqav+ktk8vdYHK2hFKO9pywvadJ+CVoTynrbXKGhFKO9py27LldVuHpwPyt13ABYBL+TrHmynraVBm+e6uRTYC7ghIo5r5yQivW5df+s3ohpyRJ3r7h9394URcWe+bSAf1KENh5wPHAasJcsb9e522lp1jb038D6yQmQH5Q0RGXP6scKlRFnt5gOBS9z9cnef5u5TgfkRsSoi1gDb5ZktWlZ6UETcNPCzu98HLC/bf/mKnQq3jYfasNB5fdixUBu2qtbw5NWzK7+jqtbwcFvbRk9LSY6oiPiMu58HfJosH9RFwMqGY9cBWwDLWj1pU5GeD5MfiYhHy/Yrq/1atX0s1IaFzuvD1r02LDRXa/j5in2aqTVcbMeWj2jzRlRpjqiIWOfupwJfzrc1/ks0jSzTRcuafQR1JPCJdk4gUgfr0oTCpURhjqjB+aDyIfHSfKg8BXg0Il5sp62VQevuhwLXRMRz7t75e34iPag/TShcilTkiBoqH9SpwCnAx4AT221r1d3j44CTgd+7+yTgc2R3k0XGlIoetVBRjqih8kFFxH3AfW2dqEHVjagLgQs7PYlIr1vX32svKxbTyxUi6I0okdppd3g8GoY1aMvmxI6Hgs7QeVHnuhd0huq5sFVzq2F45le3Qj2tSM3omlakZtTTitSM0s2I1Eyfhsci9aLhsUjNlPe0/V1rRzMUtCJASqPdguYpaEWouhGlnlak5+iaVqRm+vvbC9qixG7ufjjwUWAj4MiIiHz9F8jyRP08It7Rzjnrc59bZAT19U8oXIo0JHa7Apjl7gvz9Qa8EBELgX8BPpWv3xq4OyLmtBuwoKAVAbIbUUVLiSETu0VEiohv5+vvYn0eqLcAH3f3a929Oi9QAQWtCNDfP6FwKVGY2K3BW4HzAPIeeT7wg4F17VDQigCpZClRmtjN3bcHlkbEyxkN8174s0DbWQQVtCJA6rfCpURZYrctgd0i4j/cfYa7Tx9I9panbrqr3bYqaEWAlKxwKVKU2M3dNwNuBE539wB+RFYK5Ep3vww4hvUJ4Fo2rI98yiayj4eCztB5Uee6F3SG6gnsVQkRoLOkCPu3cUy7j3yKErsBQ/2yH9bWSQbRc1oRqBoG9xQFrQhU3nHqJQpaEdTTitRO2Q2nXqOgFQGoUdA2/cjH3V/v7teNZGNERk2bb1eMhmZLXU4GDgamj2xzREbJGLymPZqs8ldpufmyos7joaAzdF7UeSwUdK4qEP7i8+VthOoC4eVaL/uaemuee6nKoHX3twI/iYgX3L1037KCzVXbx0JBZ+i8qHP9Czo3VyB8xia/LN3eTIHwIlvTRq85xq5pjwUucPdbgN3d/cyRbZJI91l/8dJrKnvaiFg08LO73xIRZ49sk0RGQY16Wj3yEYFey91WqqWgjYgDRqgdIqOrBx/tFFFPKwLY8Cd22wFYRDYl77sR8eBQ69o5p+bTikBbL1cUJXbLnQ98FlgCLC5Z17Jh7WnL5sSOh4LO0HlR57FQ0LlqLmzV3GrocH71pk+2fIi1NzweKrHbne4+FZgfEasA3H07d585xLoNImJdqydVTysC2RtRRUuxosRus4CVDfutI8t/PHhd+dtABXRNKwLt3ogqSuz2e6DxtbJpwKoh1rX+6hbqaUWAtl+uGDKxW0SsAZa6+zR3nwI8GhHPDrHuxXbaqqAVgbZuRBUldss3nwqcAnwMOLFkXcs0PBah/Uc+RYndIuI+4L5B+75qXTsUtCL05jvGRRS0IqA3okTqRj2tSN2opxWplzbfiBoVCloRUE8rUje6phWpGQWtSN1oeCxSL+O2py2bEzseasNC5/Vh614bFqrnwlbNrYYO51cfvVHrx6inFamXcdvTitRWWU/bwlwCd58NfARYDvw8Im5v2DYR+CLwZrKJA4dHxEt5Vot7gBnA6RHx5bJzaGqeCMOarPwc4KsRcSFwurs3hvxC4EyyubebAu/M1x8JvDEi5lQFLChoRYDsjaiipUUHAw81fJ438ENE3B4RT0REH1nPuizftD/wQJ7ZsZKGxyLQVrJydz8D2GHQ6i0iYiDUB/JG/e+g4yYC0yPiVsiqeLj7LOBqd78nIm4uO2+zpS4NOIosv+u9EfF4M8eJ1EU77x5HxDmD17n7vg0fG/NGNfpz4O8HfdcKdz8F2A8oDdpmh8eLgbsi4noFrIxFwzg8vsXdt89/npwnKZ+Z9664+4HA3RHxuLtvma8buO6dBdxSdYJmSl3uQ3YB/Zi7HwF8MiJeavmPItLLhu+RzyeB4919ef4zwCfIgnki8AXgqfznq939cuBad/8WWcd4d9UJmhkeHwpcGhFfcfcvkt3OPm+oHeeWFHUeDwWdofOizvUv6FxdIPw1G0+u/I6qAuHl1rR8xHBNzYuI3wGnD1p3csPH7wxx2C6tnKOZoJ3C+iTL1wLvKtqxrGBz1faxUNAZOi/qXPeCztBcgfAHn3yhdHszBcKLzNqr9Tei6vRyRTPXtLcCe+Q/bwjcNXLNERklbaRQHS2VQRsRVwHT3X0RMBf40oi3SqTLxlQleICI+NuRbojIaOrF4CyilytEoCeHwUUUtCKA9dcnahW0IozjbIxlj23GQ0Fn6Lyoc+0LOlM9gb0qIQJ0+ghwbD/yUU8rArqmFakb9bQiNaMbUSI1M25vRInUlfWNdguap6AVgWG7EVWW2C3fPg+4newV4iMi4vvu/n6gD9gcODciSq+wlSNKhOyatmhpUVliN4BFwNw8idv38yD+ozyh2xPAYVUnUNCK0J3Ebu4+Kd++1N3/Yoj97wcOqTqBhscitPfIp9XEbnnGl4PcfRvgOne/i2xIvGLQ/qUUtCK098in3cRuEfGYu58N7Aw8RZYbqnD/wTQ8FoHhnARfldht4Bp3KvBT4EZgQb5uJ+CGqhOopxUBrG/YHtSWJXb7P+ACd78KuC3PJ4W73+Xux5ANjRdXnUBBKwLD9sinicRuew5xzJJWzqGgFUGvMYrUzrh9jbFsTux4KOgMnRd1rn1BZ6rnwlbNrYZO51evavkI9bQiNTOMN6JGnIJWBDQJXqRuNDwWqZukoBWplTF1Tevu08geFt9NVvLynIhYWX6USM3UJ2abevf47cDTEXE18Chw0Mg2SaT7rL+/cOk1zQyPfwZ8yt2vI5uFcGPRjvNK6sOOh9qw0Hl92PrXhq2uNTz7Dzer/I6qWsPl2nhOO5aGx3mZ+fOBi4ErIqKwsGhZ7deq7WOhNix0Xh+27rVhoblE47+Nx0u3N1NruMjmW7VxUI1uRFUOj919W2Br4I+Bo9x9/xFvlUi39aXipcc0MzzeC1gREWvc/XNkBaZ/NLLNEukuG6aetiyxm7vPBe4BVpIVaL8sIs4cKtlb2TmauRF1A7Ctux8CvA64vI0/i0hv6+8vXlpTlthtM2BORMwDTgW+la9/RbK3qhM0c027Gjgl/3h9C40XqY/hu6Y9GDi24fM81ueIurth/Z4R8dWGZG8nuPvJEfG1qhPo5QoR2rt73Gpit4bjpgHPw9DJ3iLiwbLzKmhFAPpafx7bbmI3spu63xv0XY3J3kqDVondRCAbHhctrSlN7JbbG7hj4MMQyd5KDWtPWzaRfTwUdIbOizrXv6Bz9QT2qoQI0FlSBF/W+ssV5TecJpZse5XCxG5kw99JwEsDQ2h334chkr2V0fBYBGCYpuZVJXbLr2HPbPh8B0MkeyujoBUB6C8rm9dbYdJbrREZLZoEL1IzPTibp4iCVgSgrz5VpRW0IlCrWT4KWhFo6+WK0aKgFQFSUtCK1It6WpGa0d1jkXpJunssUjMaHovUjG5EidSLhsciNZPG67vHi5aVz8PcuWTb01/Zo+Pzn0Ez37Fjx+cpU51f9hm2xoo3b/pkZw04eqOKHdYwa6+qfaq2V6maz7qqMjdxW3Ni27d0r5O3K8siX51Qu4ss1ej1LRFRuhmR2lHQitSMglakZhS0IjWjoBWpma48p3X3k4AngY0jYkk3zjno/DOBS8mKid0QEcd1uw15O14PnBsRfzJK5zfgKLL/F/dGVNSbHP7zTyPLVHg3sBA4JyJWdrMNY8GI97Tuvh+wWURcAcxy94Ujfc4h7A28j+xR8UHu/oZuN8DdJ5PVbJne7XM3WAzcFRHXdztgc28Hno6Iq4FHgYNGoQ21143h8SHAr/KfH8g/d1VE3BQRz+cFse8jK0PYbUcDl4zCeYGXk2IvBN7i7v+UJ83utp8Bx+QZ+GcCN45CG2qvG0G7ObAi/3mgINGoyIfJj0TEo10+71uBn+T/aIyWQ4FLI+ILwKZkNVS7Ku/dzwcuBp4Y5b+P2upG0D4FTMt/LipI1C1HkpVo6LZjyUo/3ALs7u5nVuw/EqaQFTMGuJbyt0pHhLtvC2xNVoDqKHevfutTXqUbQXs9sGv+805kRaq7zt0PBa6JiOfcvbOCOS2KiEURcUBEHEBWHfzsbp4/dyu8/HL2hsBdo9CGvYAVEbEG+FxDe6QFIx60EXEbsNrdjwaeiYgfj/Q5B3P344DPAt9x918Ao3L3djRFxFXAdHdfBMwFvjQKzbgB2NbdDwFeB1w+Cm2oPU0YEKkZvVwhUjMKWpGaUdCK1IyCVqRmFLQiNaOgFakZBa1IzShoRWrm/wFrZScyb+2aygAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOwAAADPCAYAAAD26z2MAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAY3ElEQVR4nO2deZxcVZXHv7e7s3Qn3Z1OJ53E7AQQkHU4GBbnMyAIgogRRxmGYRNxBgT3BUbRYdgyjLLIoqOAIKgz4IKCTNCRDwooyCEQWaKJQpoIgRCykK2T7q47f7zXpAj17quq113dr3O+n8/7fLrr3vfuqVd16tx3l99x3nsMw8gHdYNtgGEY5WMOaxg5whzWMHKEOaxh5AhzWMPIEeawhpEjdjiHdc4d65y7o4rzZjnnnhgIm4YCzrk659yHnXM/dc69e7DtMUoz7B3WOfd155wreukXwMRKr+O9XwasT2lrRkp5m3PuRufcPc65jzjn/tk5d1ultqS0ca5z7jPOuf9xzj1VwanHAX8FzgRW9KdNRv/RMNgGDCTOuVnACcDdRI6K937rG/23InoDbb0faAFuSarjvV/jnHsA2M17f0N83oJqjSlhQytwlPf+2Pj/Cyo4fTfgYe/9SmBlf9lk9C/D2mGB9wFnA/9C7LAxTc65K4H3AicDjwGfATYCHd77LznnzgI2AwcB/+m9/zOAc24UcDXwMHAP8BPgFOAooME59wgwDeiIr3+N9/63pYxzzjUC051zq4FvAs/G51wBHAu0AfOJegRjgX2AW4FXgDuB24Ddvfcfji+5HtjVOfdvwGXAzQn2ngT8A7AzsHd8Dw4G2p1zL8X3ayFwGFHEbQb+EZgCbPHeX+Scmxe/Pg84z3u/NOWzMPqBYdsljr+ovURf7H2cc1OKiuuATxM56UXAZOAA4AbgDufcW4EDvPc3AzcDV/Wd6L3fQvTlJ45GS+N2HgYeBJYA5wCrgd8Ae5Uwby/n3GnxdRu89+uJnPUJYH/gIaCR6EdAgVPiiDwfuNF7/xzwGnATkUP12VYA3kXkfIuJInkpewEeJf5B8t7/ichBf+69/2Ns+21EPxRTgPOBH8X36hXn3ASinssrwO8S3qMxAAznCDsPGEkUPRYBZwAXx2UbvPfeOfdr4GLv/fPOuceAp4GPEUWODXHdRcDbKmh3ItDsvV8A0WBOiTpPeu9vds7dQhSNIXL6V733vc65vr8LzrndgB4A7/0LzrkJ8TN5r/f+1eKLxj9SG7z3Rzrnjgdud84lOdPr7ZUoW0jUa9gC1BM55KbY+b/pnJsb/78AWJDwHo0BYDjf6One+yviKHku8GHnXP12daYBv3PO7QRcDxxNFEWeASSuM4boC1zMVqAp/ns80X30gANWEUX0tzvnGoiiZEl8tPNik3Nu18D7+DOwRzyKWwf8ySfv2GgEPhRf+8dEEb4jwd6SOOdagIvi+7Ypfk9LgVPj8mOBZcBxzrmdnHNjgEMD9hv9yLCMsM65jwI7O+dGe++7gBFEz4AXxc93y5xz5xB9ec8HpgPXAf9N1OV80jn3v865C4m6jZ9zzk0iet48CLgf+Jlzbg7Rc++BRE5+GVG39mzgZ0TPxqcU2dUGHALs4pz7CNH9Pw74IFEUa4ij/oHA3s65Sd77l51zlxI913YCn3XOTQOmOeeO9N4XP5sDXOWc2x14EXjCe78wfi7d3t4ZwNvjgao6oh+o+tjmeufctfHr7yPqit/pnDsRuCC26RKi7vD/AadX8TEZVeBse51h5Ifh3CU2jGGHOaxh5AhzWMPIEeawhpEjzGENI0f027TOeeedZ8PNxpBh/vz5ZS8Yn3/JuX7t+jGhKp3z58+fldmofqBf52E/dNb3E8s2r9uTxtbkzSMnfftTmdtv7iwEy3eZ0sTSFZsSy9t/m32TyuoDpwTLd57SxJ8DNoxe3ZOp/d7GcKdpzrQx/OWvG4N1xj69KpMNnX8/OVi++7gmFq9NvgcAM76+qOr29/tYuP3tWbt+DBd+4vbE8q9c/aGZVRvTzwzLhROGUSkF8tFBNIc1DKC75JLqoYc5rGEwzCKsiHyGaFNzq6peO7AmGUbt6SY8/jFUSHVYEXkH0K6qXxORC0Rkrqo+UgPbDKNm9Fa5pj4pmInIrkR7hjcBd6nqEhE5iWgzyEjgGFXdHNe9k2hDxl2qeub2bRRTzjzsMUSboSHakXJMZW/JMIY+3fjEI4miYHYr0CYic4uKrwauBK4l2u0EsFBVjyLarrhLfI0DgG+o6uQ0Z4XyusQTgDXx311E6gwl2bxuz8SLbNkUHhnfvbUpWF4OjVPC3ZopbaOC5c17jM9sw/gp4fcxZXzYhpHN2bpmhZHh6cfJE0anXmN0IdsATNO48D2YNiZ8DwAm7FfZ1ExWuqsLsKWC2SMi0gjMUdUNACIyW0QaVHWxiDgiVZK+Oc7DgHNF5D7gLFUNzneV47CvsG3zczPwalLF0DxrWvnide8qw5QwzSvSv+zBedhnVme2YXVL+pdxMOdhgfR52MVrguVpdL4t/R6kzsM+/lLV7Y8/uHJn76UqYb6kYNZGJOHTRw+REskKor3SnwUeAB5R1ctF5ArgP4DzgC+HGiynS3wPkVAXwB5Av6n8GcZQodu7xCNAUjB7FSjuyjQBawFU9Trg4xQJG6hqD/AFYHaanakOq6oPAV0icjqwVlV/k3aOYeSNXlziEWD7YHaviLSq6hagU0SaRGQ0sLxvgCnmWSL9MOIuMkQO/2CanWVN66jqxem1DCO/dPvK98Go6kMiclhfMIuPbwInEkXMzxMJ2X1aRJqJZIN+SCQ79O34Mg+KyOPA40SqnUFs4YRhAL1VblwrEcxOjF9/im0DS30cVuL8QyppzxzWMKguwg4G/eqwZ170ycSyXSc1seTlIxLLx1Y5rl5M6/fD6zmaZCqt+kJieW/DiMw2jH8gPMrbvGc7459KHGine2bFaX/eQNPTrwTLG7vaaQ60D9DzwouZbJh22bPB8naZyrTA5wDQ+aWDq25/v60vV3xO747osIaRV7rZXrJ6aGIOaxhAtzeHNYzcUO2gU60xhzUMoNvnwxXyYaVhDDC94RVNQwZzWMPAIqxh5Ap7hjWMHGGjxIaRI8xhDSNH2Eonw8gRFmENI0fYoJNh5IhqI2xW1cRS9ULt5eNnxTAGmIKvSzyS6A/VxIR6iZjDGgZRhE06ApSUAC5WTYzlYkqqJibVCzVoDmsYVO2wlaomwjbVxANS6pWkX59hGzYnb0Kv2+KD5a0/XJi5/cJBewfL/Zxm/Ij2xPL6pdk2bgMUVoclQgvrRwXrNGzZms2A0SkSo85FR6jKyJGZTKhraQmXt7ZSPzH8Pmd/Z1n1BpzUWPEpVU7rVKWaKCKvED3LfjqpXhIWYQ2DqiNsJtXEMuq9CRslNgygENitkxTV+kk18Q310uw0hzUMwtM6oYeMflBNLFUvEXNYwyAcYYcS5rCGQTjCVj6ENXAEHTbud98E7A8sUNWza2KVYdSYnkI+1hKnjRIfCJwG7AkcHueyNIxhRwGXeAwlghFWVX/Z97eIPAUEcwDuHMiNmpYXdcz+U4Ll5VCY0xwsnzQ53Lmpb8wm4g3gu8I2dOzUFix3o9PztwbPHxkWQ++YEbYPoLA225OSawrnh+2Y1Zp+jYYsES+cyrIU3TmJsGV9MnHX+HlVXR6qF8p7mlbe+tiKckwJUqgbl1pn2V/WJ5Y1LA2r5peD3xjOvQrQGch96saMydS+S1s4AXSmKP/3rsx2H9IWTgB0/iHcRtoPT4hx+1T+1DncBp1OJiXRrGHkmZ6cbGBPtVJE5gF3qup6EZlUA5sMo+ZUs1tnMEgbJT4b+BzwqoiMBK4iGjU2jGFFXiJs2qDT9cD1NbLFMAaNnsIwcFjD2FEYboNOhjGsGRZd4kppuy85ke/YvSbS9mTyUH5vd8Z9oEDDM8uC5fUjJtHwTOXJfiuhrjVlL+iYpmAdvzXbfehZ/tdgee8kT8/ycDLluub0udoQ3btNC9swbSzdW8PTT5s7suzJXVXxGRZhDSNH2DOsYeSIaiNshaqJJwKfBFqAk1VV47rXAB8EnlDVd4fay8fPimEMML2+LvFIohLVxFh8bZOqzgW+ClwYX2MqkZri5DRnBXNYwwCgt1CXeAQoWzURqFfVn8Z1HwX61uK+E7hARO4WkQlpdprDGgZRlzjpCFCNaiLAEcAVAHF0ngP8qu+1EOawhkHVEbZi1UQR2RnoVNVn+gpV1avqlUQZAYLYoJNhAD5ZgTfEPcDRwO28UTVxnYh0ikgTUCBWQ4zX4u+jqj8SkbGAJ3qu9fHS30fTGjSHNQzSdIkLJV+tUDWxHbgX6BGR8wEHCHC7iGwEHonPDWIOaxhUP61ToWriviUu8cFK2jOHNQygULCVToaRG1IGl4YM5rCGQdWDTjXHHNYwgIJFWMPIDzkJsOawhgHgbdDJMPKD3xH3w65875zEsnEdTayclbxxu+Ou7O0XXl0dLPddXRQ2JOsGZ924DVB4LVn3GKCwaUywTmFzV6b2G6ZMDpbXj2+jYUpvsI7vCZen0T02/LXqGV2fWqerrbbPlDatYxg5wrrEhpEncjLqZA5rGFiENYxcsUMOOhlGbsmJw5Y9FCciu4nIzwfSGMMYNHzgGEKUm25yFHAkkC0XomEMVYbZM+zpwA3A8aFKu3YkJ/J9S2tYOLp1r35Iprw+nAy5Y+fxwXLXWHle0TfZ0NMTtiEloXNWIfH6ceFkyR0z03O3+t7SG7bLZcuM8O/65InpSau72jMkWA5/BCXxVb7lrDKnItIBnEOULP0JVf1tqL1UhxWRI4AHVHWTiATrLlkZTugcKu8IZAUol7SFEwDLHnsxsaw/Fk6U43ChhM6ZF05M6k5v/8mwMn7WhRObRqZ3xJ59Ppz4ekNv9cMrs8K/y6Wp4hm2SOb0ayJygYjMVdVH4uKriTandwM/EJEPEMucisgZRDKn7wEuBS6PHfouETlOVRM74uU8w54JXCci9wP7isgXK35nhjHEcYXkI0B/yJweCSwtuuasUIOpP2OqekLf3yJyv6peknaOYeSO6kaJK5U57XPS12VOgRFFEbXvGs8lNWjTOoYBSTprafSHzOmGonrF1yhJRSusVfXQSuobRm6oblrnHmDv+O9imdMtQKeINInIaBJkTkVkDHB/7MQAo1R1SajBfGyzN4wBxhVc4pGEqj4EdJWQOYVtMqef4o0yp+eLiAK/JhpB/gpwhoh8Iv47iHWJDQOqXiDRDzKnLwLnl9tevzrsyA3J77qhxQfLu3ebmrn9+ofSp3VCFNaH97KWQ9o8rO/qorAxPKWRicaUOc6RI9PrrHstXJ5CU+e6YPmoRpdaZ9Ok9kw2VIobYiuakrAIaxgw7FY6GcbwxiKsYeSHlAUSQwZzWMMAi7CGkSdC0zdDCXNYw8C6xIaRL6xLbBj5wSKsYeQJi7CGkR9spZNh5AlzWMPID/YMaxg5whzWMPJElV3iClUT24AvAY+r6m1xvTrgYWAGcL2q/nuoPdvAbhhUJ8JWpJp4K9AmInOLiq8GrgSuBeYDqOoa4E+8MVAeD5yiqpPTnBX6OcLWb03+marrCZePeDG8P7IsUmRKXWNjUMrUb9mS2YT6aW8Jl08aT8OsZM1dPyLjR5ImUep9dITorkLYt5itKVKr3T2pddrv+EP17X8snCO3JNVF2FKqiY8UqyYCiMhsEWlQ1R5gex3cQ4BrROR7wOdVNdg5twhrGFQtc1qpauKbUNVPAXOAacBpaXaawxoGVCvCVrFqYilUdRPwcUpLyLwBc1jDoOoIW5FqYqkLiEjfNqF24FdpdprDGgbRSqekI4lKVBMBYllTAf4mduYmYJGIXA7sW5QZIBGb1jEMqFZIvCLVRFXdSJT4qpi9qYBy00064FSi+aZFqvpCJY0YxlBnuK0lng98V1WfHkhjDGOwyIvDpj7DishBwFzgnSJymYiMHHizDKPGFALHEKKcCDsPuElVvysi/0XUB7+iVMU505Lzgk5pDyd0bkwR4C6L10YEiztmh5Md++703Kpp1I0LJ0zumBFe3OHrMyQyBlxKMuaOGWNTr+HXZ/xNbg2/x47pZdgwubbjoXmJsOU47Gi2TQLfDbw/qeJf/hpWtA+VNy9ek1hWNqvSr9G5aGViWX+sdKrrSP/h6Xw6OUNB1pVOroxkzJ3PhDMk+Fczfhbl3IOUz9uvSP6c0hh/cOUrnfKy+L+cn7EHgf3iv0cQJaM1jOFFdQsnak6qw6rqHcAYETkBmAncOOBWGUaNqXLhRM0pq/+lqp8daEMMYzAZao6ZhC2cMAwYcl3fJMxhDQNwhXx4rDmsYTC8pnXKpmXhisSypi3jaQlMZ/Qsez5z+64h/Hb85uZg0ub6Cf2QRDhtc7hPqbMm20Z+nzIP6zeOxq8NJ2wu7DI9kw3uj8vCNkyuS522eem0fTJY8HLFZ9gzrGHkiR0xwhpGXrEIaxg5otpBp35QTTyYSNepDviOqgafFWwDu2FQ3Qb2flJNvAz4KvAD4MI0O81hDQNwvclHgFKqiRSrJsZyMbNFpM9JX1dNjKNwj6p6VX0e+Ns0O81hDQOqXUucVTWx+Py+84LYM6xhUPUzbFbVxOLzAVK3i1mENQyqe4Ylo2qiqi4ldmwR2Qm4P81Oc1jDoLrdOllVE+N6F4vI54F/Ar6YZqd1iQ2D6qd1sqomqup9wH3ltmcOaxhgK50MI0+43nx4rDmsYYBFWMPIE7Yf1jByxA65H/bZU6cllo1qaeJZaUosn3Fh9v2wdU3J1wdwo0YF6/j2cZlt6J4QtqGnYyxbN7rE8vrXts/3WxmuN7yWzk9oxU9Pbh+gbunyTDYsuXCvYPnIpiaWziuZLvV12hdl8KBkeexELMIaRo6wQSfDyBP58FdzWMMA6xIbRr5I0+IaIpjDGgbD6Bk2XqR8PrCQKO3kpaoalt0zjLyRD38ta7fOUcAqVf0JsBw4fGBNMoza4wqFxGMoUU6X+PfAhSLyc6JNuvcmVdy9JXkOcnpTOD/sRJlahilh0uZhO3ZKmWedEM4fWw69LeH3OWlSY7C8bnPGPLlp+WGnpk9SunHpKStDjEj5HGaMCt8jgObJ4bni/mbYdIlV9QURuRr4FnCrqm5Kqrv4tcSi1PKN+kKaKanUt4STKQN0PhEQmZ6efXtw2sIJgOee25BYNtALJwA6l4TFyt2yyoW4i1lyXLog+9Obwt+V9peqd9gZU6o4qcpBpwpVE9/0Wlz3TuDA+LUzQ+2lfkNFZDowFTgaOFVE/q6qd2YYQ5len3wkUKlqYqnXROQA4BuqOjnNWaG8Z9j9gTWx7MVVbEvubBjDBud94hGgEtXE5hKvNQCHATeIyC1FKhSJlOOwC4DpInIM8FbgljLOMYx8USgkH8lUoprYUuK1iap6OTAbWAWcl2ZmOc+wXUTaNBCJThnG8KO6Z9hKVBM3lHhtLYCq9ojIF4DvpDVoImyGQTRKnHQEqEQ1cV2J1zaLSN/oWjPwYJqdttLJMCB1OqwUqvqQiBxWQjXxRLapJm4hVk1MeO1BEXkceBy4Ia1Nc1jDgJQucfIUU4WqiaVeO6QSM/vVYWd/Lzmh86Td2+havCax/KWzDsrc/oRF4bm93pnN9PQkD8StOKSKnc/btzEyXN7SOooXpibPlY5ZMTqxrBwm/viZYLlrK+CeD8+zrn7vHplsmPGL8OKPiTMLzOgM17nv5tRgk8i/XvWByk8KDi7VV21Lf2MR1jAAbHudYeSIQmiF2NBxk6FjiWEMJhZhDSNHDLFdOUmYwxoGQBmbJoYC5rCGASYRYxi5ooqFE4OBOaxhAN6bwxpGfrAIaxg5wkaJDSM/eBslNowcYV1iw8gRNuhkGPnBusSGkSP8jriWeN/jQ4oz62jbPVS+MrsB70irsIEJQc3aZL3gsulOKV8VKdklMiFj+x9Nu0AvbXPT6rySzYZUTfg1dLSFa1S1p7V6Ovf/3OyZofKaWZKC8zlZkmUYhomwGUauMIc1jBxhDmsYOcIc1jByhDmsYeSImszDJqXkqxVxIqKbiBJ7LVDVs2ttQ2zHbsDXVPU9g9S+A04l+iwWqfZDjs/K2m8CzgcWAnOBS1X1tfBZRjEDHmFTUvLVigOB04A9gcPjFH81RURGAUcC2cWPq2c+8Kiq3lNrZ405Clilqj8BlgOHD4INuaYWXeKSKflqiar+UlU3xsmonwJeqrUNwOmUkYphoBCRg4ii2jtF5DIRSZE8HxB+D5whIjsT5ZK5dxBsyDW1cNiklHw1J+4aP6+qy2vc7hHAA6Hs9TVgHnCTql4DjAfOqbUBcVS/GvgW8PIg349cUguHTUrJNxicDHx5ENo9E7hORO4H9hWRLw6CDaPZlp/0bqLHg5oiItOJFi4eDZwqIn9XaxvyTi0cdvuUfAtq0OabEJF5wJ2qul5EJtWybVU9QVUPVdVDgSdU9ZJath/zILBf/PcI4NFBsGF/YE2cjvGqInuMMhlwh1XVh4CuvpR8qvqbgW5ze0TkbOBK4Gci8gdgUEZpBxNVvQMYIyInADOBGwfBjAXAdBE5hmgPxC2DYEOuscX/hpEjbOGEYeQIc1jDyBHmsIaRI8xhDSNHmMMaRo4whzWMHGEOaxg5whzWMHLE/wMwT1YhihpipAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Set model and likelihood into evaluation mode\n",
    "model.eval()\n",
    "likelihood.eval()\n",
    "\n",
    "# Generate nxn grid of test points spaced on a grid of size 1/(n-1) in [0,1]x[0,1]\n",
    "n = 10\n",
    "test_x = torch.zeros(int(pow(n, 2)), 2)\n",
    "for i in range(n):\n",
    "    for j in range(n):\n",
    "        test_x[i * n + j][0] = float(i) / (n-1)\n",
    "        test_x[i * n + j][1] = float(j) / (n-1)\n",
    "\n",
    "with torch.no_grad(), gpytorch.settings.fast_pred_var():\n",
    "    observed_pred = likelihood(model(test_x))\n",
    "    pred_labels = observed_pred.mean.view(n, n)\n",
    "\n",
    "# Calc abosolute error\n",
    "test_y_actual = torch.sin(((test_x[:, 0] + test_x[:, 1]) * (2 * math.pi))).view(n, n)\n",
    "delta_y = torch.abs(pred_labels - test_y_actual).detach().numpy()\n",
    "\n",
    "# Define a plotting function\n",
    "def ax_plot(f, ax, y_labels, title):\n",
    "    if smoke_test: return  # this is for running the notebook in our testing framework\n",
    "    im = ax.imshow(y_labels)\n",
    "    ax.set_title(title)\n",
    "    f.colorbar(im)\n",
    "\n",
    "# Plot our predictive means\n",
    "f, observed_ax = plt.subplots(1, 1, figsize=(4, 3))\n",
    "ax_plot(f, observed_ax, pred_labels, 'Predicted Values (Likelihood)')\n",
    "\n",
    "# Plot the true values\n",
    "f, observed_ax2 = plt.subplots(1, 1, figsize=(4, 3))\n",
    "ax_plot(f, observed_ax2, test_y_actual, 'Actual Values (Likelihood)')\n",
    "\n",
    "# Plot the absolute errors\n",
    "f, observed_ax3 = plt.subplots(1, 1, figsize=(4, 3))\n",
    "ax_plot(f, observed_ax3, delta_y, 'Absolute Error Surface')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## KISS-GP for higher dimensional data w/ Additive Structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The method above won't scale to data with much more than ~4 dimensions, since the cost of creating the grid grows exponentially in the amount of data. Therefore, we'll need to make some additional approximations.\n",
    "\n",
    "If the function you are modeling has additive structure across its dimensions, then SKI can be one of the most efficient methods for your problem.\n",
    "\n",
    "To set this up, we'll wrap the `GridInterpolationKernel` used in the previous two models with one additional kernel: the `AdditiveStructureKernel`. The model will look something like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "class GPRegressionModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super(GPRegressionModel, self).__init__(train_x, train_y, likelihood)\n",
    "        \n",
    "        # SKI requires a grid size hyperparameter. This util can help with that\n",
    "        # We're setting Kronecker structure to False because we're using an additive structure decomposition\n",
    "        grid_size = gpytorch.utils.grid.choose_grid_size(train_x, kronecker_structure=False)\n",
    "        \n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = gpytorch.kernels.AdditiveStructureKernel(\n",
    "            gpytorch.kernels.ScaleKernel(\n",
    "                gpytorch.kernels.GridInterpolationKernel(\n",
    "                    gpytorch.kernels.RBFKernel(), grid_size=128, num_dims=1\n",
    "                )\n",
    "            ), num_dims=2\n",
    "        )\n",
    "\n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n",
    "\n",
    "    \n",
    "likelihood = gpytorch.likelihoods.GaussianLikelihood()\n",
    "model = GPRegressionModel(train_x, train_y, likelihood)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Essentially, the `AdditiveStructureKernel` makes the base kernel (in this case, the `GridInterpolationKernel` wrapping the `RBFKernel`) to act as 1D kernels on each data dimension. The final kernel matrix will be a sum of these 1D kernel matrices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scaling to more dimensions (without additive structure)\n",
    "\n",
    "If you can't exploit additive structure, then try one of these other two methods:\n",
    "\n",
    "1. SKIP - or [Scalable Kernel Interpolation for Products](./Scalable_Kernel_Interpolation_for_Products_CUDA.ipynb)\n",
    "2. KISS-GP with [Deep Kernel Learning](../06_PyTorch_NN_Integration_DKL/KISSGP_Deep_Kernel_Regression_CUDA.ipynb)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
